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ABSTRACT 

The Jeans equations relate the second-order velocity moments to the density and po- 
tential of a stellar system. For general three-dimensional stellar systems, there are three 
equations and six independent moments. By assuming that the potential is triaxial 
and of separable Stackel form, the mixed moments vanish in confocal ellipsoidal coor- 
dinates. Consequently, the three Jeans equations and three remaining non- vanishing 
moments form a closed system of three highly-symmetric coupled hrst-order partial 
differential equations in three variables. These equations were first derived by Lynden- 
Bell, over 40 years ago, but have resisted solution by standard methods. We present 
the general solution here. 

We consider the two-dimensional limiting cases first. We solve their Jeans equa- 
tions by a new method which superposes singular solutions. The singular solutions, 
which are new, are standard Riemann-Green functions. The resulting solutions of the 
Jeans equations give the second moments throughout the system in terms of pre- 
scribed boundary values of certain second moments. The two-dimensional solutions 
are applied to non-axisymmetric discs, oblate and prolate spheroids, and also to the 
scale- free triaxial limit. There are restrictions on the boundary conditions which we 
discuss in detail. We then extend the method of singular solutions to the triaxial case, 
and obtain a full solution, again in terms of prescribed boundary values of second 
moments. There are restrictions on these boundary values as well, but the boundary 
conditions can all be specified in a single plane. The general solution can be expressed 
in terms of complete (hyper)elliptic integrals which can be evaluated in a straightfor- 
ward way, and provides the full set of second moments which can support a triaxial 
density distribution in a separable triaxial potential. 

Key words: celestial mechanics, stellar dynamics - galaxies: elliptical and lenticular, 
cD - galaxies: kinematics and dynamics - galaxies: structure 



1 INTRODUCTION 

Much has been learned about the mass distribution and in- 
ternal dynamics of galaxies by modeling their observed kine- 
matics with solutions of the Jeans equations (e.g., Binney & 
Tremaine 1987). These are obtained by taking velocity mo- 
ments of the collisionless Boltzmann equation for the phase- 
space distribution function /, and connect the second mo- 
ments (or the velocity dispersions, if the mean streaming 
motion is known) directly to the density and the gravita- 
tional potential of the galaxy, without the need to know 
/. In nearly all cases there are fewer Jeans equations than 
velocity moments, so that additional assumptions have to 
be made about the degree of anisotropy. Furthermore, the 
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resulting second moments may not correspond to a physi- 
cal distribution function / > 0. These significant drawbacks 
have not prevented wide application of the Jeans approach 
to the kinematics of galaxies, even though the results need 
to be interpreted with care. Fortunately, efficient analytic 
and numerical methods have been developed in the past 
decade to calculate the full range of distribution functions / 
that correspond to spherical or axisymmetric galaxies, and 
to fit them directly to kinematic measurements (e.g., Ger- 
hard 1993; Qian et al. 1995; Rix et al. 1997; van der Marel et 
al. 1998). This has provided, for example, accurate intrinsic 
shapes, mass-to-light ratios, and central black hole masses 
(e.g., Verolme et al. 2002; Gebhardt et al. 2003). 

Many galaxy components are not spherical or axisym- 
metric, but have triaxial shapes (Binney 1976, 1978). These 
include early-type bulges, bars, and giant elliptical galax- 
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ies. In this geometry, there are three Jeans equations, but 
little use has been made of them, as they contain six inde- 
pendent second moments, three of which have to be chosen 
ad-hoc (see, e.g., Evans, Carollo & de Zeeuw 2000). At the 
same time, not much is known about the range of physical 
solutions, as very few distribution functions have been com- 
puted, and even fewer have been compared with kinematic 
data (but see Zhao 1996). 

An exception is provided by the special set of triaxial 
mass models that have a gravitational potential of Stackel 
form. In these systems, the Hamilton-Jacobi equation sep- 
arates in orthogonal curvilinear coordinates (Stackel 1891), 
so that all orbits have three exact integrals of motion, which 
are quadratic in the velocities. The associated mass distri- 
butions can have arbitrary central axis ratios and a large 
range of density profiles, but they all have cores rather than 
central density cusps, which implies that they do not provide 
perfect fits to galaxies (de Zeeuw, Peletier & Franx 1986). 
Even so, they capture much of the rich internal dynamics 
of large elliptical galaxies (de Zeeuw 1985a, hereafter Z85; 
Statler 1987, 1991; Arnold, de Zeeuw & Hunter 1994). Nu- 
merical and analytic distribution functions have been con- 
structed for these models (e.g., Bishop 1986; Statler 1987; 
Dejonghe & de Zeeuw 1988; Hunter & de Zeeuw 1992, here- 
after HZ92; Mathieu & Dejonghe 1999), and their projected 
properties have been used to provide constraints on the 
intrinsic shapes of individual galaxies (e.g., Statler 1994a, 
b; Statler & Fry 1994; Statler, Dejonghe & Smecker-Hane 
1999; Bak & Statler 2000; Statler 2001). 

The Jeans equations for triaxial Stackel systems have 
received little attention. This is remarkable, as Eddington 
(1915) already knew that the velocity ellipsoid in these mod- 
els is everywhere aligned with the confocal ellipsoidal coor- 
dinate system in which the motion separates. This means 
that there are only three, and not six, non-vanishing second- 
order velocity moments in these coordinates, so that the 
Jeans equations form a closed system. However, Eddington, 
and later Chandrasekhar (1939, 1940), did not study the 
velocity moments, but instead assumed a form for the dis- 
tribution function, and then determined which potentials 
are consistent with it. Lynden-Bell (1960) was the first to 
derive the explicit form of the Jeans equations for the triax- 
ial Stackel models. He showed that they constitute a highly 
symmetric set of three first-order partial differential equa- 
tions (PDEs) for three unknowns, each of which is a function 
of the three confocal ellipsoidal coordinates, but he did not 
derive solutions. When it was realized that the orbital struc- 
ture in the triaxial Stackel models is very similar to that in 
the early numerical models for triaxial galaxies with cores 
(Schwarzschild 1979; Z85), interest in the second moments 
increased, and the Jeans equations were solved for a number 
of special cases. These include the axisymmetric limits and 
elliptic discs (Dejonghe & de Zeeuw 1988; Evans & Lynden- 
Bell 1989, hereafter EL89), triaxial galaxies with only thin 
tube orbits (HZ92), and, most recently, the scale-free limit 
(Evans et al. 2000). In all these cases the equations simplify 
to a two-dimensional problem, which can be solved with 
standard techniques after recasting two first-order equations 
into a single second-order equation in one dependent vari- 
able. However, these techniques do not carry over to a single 
third-order equation in one dependent variable, which is the 



best that one could expect to have in the general case. As a 
result, the general case has remained unsolved. 

In this paper, we first present an alternative solution 
method for the two-dimensional limiting cases which does 
not use the standard approach, but instead uses superposi- 
tions of singular solutions. We show that this approach can 
be extended to three dimensions, and provides the general 
solution for the triaxial case in closed form, which we give 
explicitly. We will apply our solutions in a follow-up paper, 
and will use them together with the mean streaming mo- 
tions (Statler 1994a) to study the properties of the observed 
velocity and dispersion fields of triaxial galaxies. 

In we define our notation and derive the Jeans 
equations for the triaxial Stackel models in confocal ellip- 
soidal coordinates, together with the continuity conditions. 
We summarise the limiting cases, and show that the Jeans 
equations for all the cases with two degrees of freedom cor- 
respond to the same two-dimensional problem. We solve this 
problem in SjJ] first by employing a standard approach with 
a Riemann-Green function, and then via the singular so- 
lution superposition method. We also discuss the choice of 
boundary conditions in detail. We relate our solution to that 
derived by EL89 in Appendix^ and explain why it is dif- 
ferent. In we extend the singular solution approach to 
the three-dimensional problem, and derive the general solu- 
tion of the Jeans equations for the triaxial case. It contains 
complete (hyper)elliptic integrals, which we express as single 
quadratures that can be numerically evaluated in a straight- 
forward way. We summarise our conclusions in J^j 



2 THE JEANS EQUATIONS FOR SEPARABLE 
MODELS 

We first summarise the essential properties of the triaxial 
Stackel models in confocal ellipsoidal coordinates. Further 
details can be found in Z85. We show that for these models 
the mixed second-order velocity moments vanish, so that the 
Jeans equations form a closed system. We derive the Jeans 
equations and find the corresponding continuity conditions 
for the general case of a triaxial galaxy. We then give an 
overview of the limiting cases and show that solving the 
Jeans equations for the various cases with two degrees of 
freedom reduces to an equivalent two-dimensional problem. 

2.1 Triaxial Stackel models 

We define confocal ellipsoidal coordinates (A, fi, v) as the 
three roots for r of 

2 2 2 

1 - 1 = 1, (2.1) 

T + QT + /3T+7 

with (x,y, z) the usual Cartesian coordinates, and with con- 
stants a, /3 and 7 such that —7 < v < —f3 < /1 < —a < A. 
For each point (x,y,z), there is a unique set (A,/x, v), 
but a given combination (A, [i, v) generally corresponds to 
eight different points (±x, ±y, ±z). We assume all three- 
dimensional Stackel models in this paper to be likewise eight- 
fold symmetric. 

Surfaces of constant A are ellipsoids, and surfaces of 
constant p and v are hyperboloids of one and two sheets, 
respectively (Fig. 0. The confocal ellipsoidal coordinates 
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are approximately Cartesian near the origin and become a 
conical coordinate system at large radii with a system of 
spheres together with elliptic and hyperbolic cones (Fig-EJ- 
At each point, the three coordinate surfaces are perpendic- 
ular to each other. Therefore, the line element is of the form 
ds 2 — P 2 d\ 2 + Q 2 dp 2 + R 2 du 2 , with the metric coefficients 



Q 



4(A + a)(A + /3)(A + 7 )' 
{p-v)(p-\) 

4(/i + a)(p + j8)0i + 7)' 
(u-\)(v- p) 



(2.2) 



We restrict attention to models with a gravitational poten- 
tial Vs(X,/J,,i>) of Stackel form (Weinacht 1924) 



V S 



F(X) 



F[u) 



(2.3) 



(X-fj,)(X-v) (ji-v)((i-X) {y-X){y-nY 

where F(r) is an arbitrary smooth function. 

Adding any linear function of r to F(t) changes Vs by 
at most a constant, and hence has no effect on the dynamics. 
Following Z85, we use this freedom to write 



F(t) = (t + o)(t + 7 )G(t), 



(2.4) 



where G(r) is smooth, ft equals the potential along the in- 
termediate axis. This choice will simplify the analysis of the 
large radii behaviour of the various limiting cases. 1 

The density ps that corresponds to Vs can be found 
from Poisson's equation or by application of Kuzmin's 
(1973) formula (see de Zeeuw 1985b). This formula shows 
that, once we have chosen the central axis ratios and the 
density along the short axis, the mass model is fixed ev- 
erywhere by the requirement of separability. For centrally 
concentrated mass models, Vs has the x-axis as long axis 
and the z-axis as short axis. In most cases this is also true 
for the associated density (de Zeeuw et al. 1986). 

2.2 Velocity moments 

A stellar system is completely described by its distribution 
function (DF), which in general is a time-dependent func- 
tion / of the six phase-space coordinates (x, v). Assuming 
the system to be in equilibrium (df/dt = 0) and in steady- 
state (df/dt — 0), the DF is independent of time t and 
satisfies the (stationary) collisionless Boltzmann equation 
(CBE). Integration of the DF over all velocities yields the 
zeroth-order velocity moment, which is the density p of the 
stellar system. The first- and second-order velocity moments 
are defined as 



<««>(*) = - 



«i/(x, v) d v, 



ViVjf(x,v) d v, 



(2.5) 



where i,j — 1,2,3. The streaming motions (vi) together 
with the symmetric second-order velocity moments (viVj) 
provide the velocity dispersions a 2 j = (i>ii>j) — (vi){vj). 

1 Other, equivalent, choices include F(t) = — (r + ct)(r + 7)G(r) 
by HZ92, and F(t) = (t + a)(r + P)V{t) by de Zeeuw et al. 
(1986), with U(t) the potential along the short axis. 




Figure 1. Confocal ellipsoidal coordinates. Surfaces of constant 
A are ellipsoids, surfaces of constant p are hyperboloids of one 
sheet and surfaces of constant v are hyperboloids of two sheets. 



The continuity equation that results from integrating 
the CBE over all velocities, relates the streaming motion 
to the density p of the system. Integrating the CBE over 
all velocities after multiplication by each of the three veloc- 
ity components, provides the Jeans equations, which relate 
the second-order velocity moments to p and V, the poten- 
tial of the system. Therefore, if the density and potential 
are known, we in general have one continuity equation with 
three unknown first-order velocity moments and three Jeans 
equations with six unknown second-order velocity moments. 

The potential 12.31 is the most general form for which 
the Hamilton-Jacobi equation separates (Stackel 1890; 
Lynden-Bell 1962b; Goldstein 1980). All orbits have three 
exact isolating integrals of motion, which are quadratic in 
the velocities (e.g., Z85). It follows that there are no irreg- 
ular orbits, so that Jeans' (1915) theorem is strictly valid 
(Lynden-Bell 1962a; Binney 1982) and the DF is a func- 
tion of the three integrals. The orbital motion is a com- 
bination of three independent one-dimensional motions — 
either an oscillation or a rotation — in each of the three 
ellipsoidal coordinates. Different combinations of rotations 
and oscillations result in four families of orbits in triaxial 
Stackel models (Kuzmin 1973; Z85): inner (I) and outer (O) 
long-axis tubes, short (S) axis tubes and box orbits. Stars 
on box orbits carry out an oscillation in all three coordi- 
nates, so that they provide no net contribution to the mean 
streaming. Stars on I- and O-tubes carry out a rotation in 
v and those on S-tubes a rotation in p, and oscillations in 
the other two coordinates. The fractions of clockwise and 
counterclockwise stars on these orbits may be unequal. This 
means that each of the tube families can have at most one 
nonzero first-order velocity moment, related to p by the con- 
tinuity equation. Statler (1994a) used this property to con- 
struct velocity fields for triaxial Stackel models. It is not 
difficult to show by similar arguments (e.g., HZ92) that all 
mixed second-order velocity moments also vanish 



{vW/j) = (VpVv) = {v v v\) = 0. 



(2.6) 
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Eddington (1915) already knew that in a potential ol the 
form the axes ol the velocity ellipsoid at any given 

point are perpendicular to the coordinate surfaces, so that 
the mixed second-order velocity moments are zero. We are 
left with three second-order velocity moments, (v x ), {v 2 ^) 
and {vl), related by three Jeans equations. 

2.3 The Jeans equations 

The Jeans equations for triaxial Stackel models in confo- 
cal ellipsoidal coordinates were first derived by Lynden-Bell 
(1960). We give an alternative derivation here, using the 
Hamilton equations. 

We first write the DF as a function of (A, p, u) and the 
conjugate momenta 



2 d\ 2 dp, 2 dv 

Px = P—, Pv = Q—, Vv = R 



dt 1 r " "* dv r " " dt 1 (2 ' 7) 

with the metric coefficients P, Q and R given in 12. 21 . In 
these phase-space coordinates the steady-state CBE reads 



dr df_ dp^df_ _ Q 
dt dr dt dp T 



(2.8) 



(2.9) 



where we have used the summation convention with respect 
to t = A, p, v. The Hamilton equations are 

dr_dH_ dp T _ dH 

dt dp T ' dt dr ' 

with the Hamiltonian defined as 

H =3* + w + & +vix >^- (2 - l0) 

The first Hamilton equation in 12.91 defines the momenta 
tfF?l and gives no new information. The second gives 



dpx p|ap ridQ ,rt9R_dv 

dt P 3 dX Q 3 dX R 3 dX dX ' 



(2.11) 



and similar for p M and p v by replacing the derivatives with 
respect to A by derivatives to p and v, respectively. 

We assume the potential to be of the form Vs defined in 
lEHt . and we substitute <l2~7t and ||2~TO in the CBE J2"5l . 
We multiply this equation by p x and integrate over all mo- 
menta. The mixed second moments vanish 12.61 . so that we 
are left with 

Hfpl) dP , (fpl) dQ + {fpl) OR 



P 3 dX 



Q 3 dX R 3 dX 
1 d 



/ , ()A (M>-</)^=0, (2.12) 



where we have defined the moments 



fd"p = PQRp, 



(2.13) 



(fpl) = J plfd A p = P s QRT xx , 

with the diagonal components of the stress tensor 

T rT (X,fi,^) = p{v 2 T ), T = X,tx,v. (2.14) 

The moments (fpfy and (fpt) follow from {fp\) by cyclic 
permutation A — *• fi — > v — > A, for which P^Q — > R—> P. We 
substitute the definitions 12.131 in eq. lETTil and carry out 
the partial differentiation in the fourth term. The first term 




Figure 2. Special surfaces inside (X = —ct) and outside (fi = — ct) 
the focal ellipse in the plane x = 0, and inside (fi = —(3) and 
outside (v = — j3) the two branches of the focal hyperbola in the 
plane y = and the plane 2 = (v = —-f). 

in 12.121 then cancels, and, after rearranging the remaining 
terms and dividing by PQR, we obtain 

dT xx Txx-T^ dQ T xx -T vv dR _ 8V S . 

~dT + q 9A + r ~dx~~ p ^x- (2 - 15) 

Substituting the metric coefficients 12.21 and carrying out 
the partial differentiations results in the Jeans equations 



dT x 



dX 



+ 



T\ \ — T u 



+ 



T xx — T„ 



2(A-m) 2{X-v) 



+ 



T 



drt 2(n-v) 2(/i-A) 



OTi/is T uu T\\ Tuu 



= - P % (2.16a) 
= -p^, (2.16b) 



-p— , (2.16c) 



dv 2(u-X) 2{v-p) 

where the equations for /j, and v follow from the one for A 
by cyclic permutation. These equations are identical to those 
derived by Lynden-Bell (1960). 

In self-consistent models, the density p must equal ps, 
with ps related to the potential Vs 12.31 by Poisson's equa- 
tion. The Jeans equations, however, do not require self- 
consistency. Hence, we make no assumptions on the form 
of the density other than that it is triaxial, i.e., a function 
of (A, p., v), and that it tends to zero at infinity. The result- 
ing solutions for the stresses T TT do not all correspond to 
physical distribution functions / > 0. The requirement that 
the T TT are non-negative removes many (but not all) of the 
unphysical solutions. 

2.4 Continuity conditions 

We saw in H2.2I that the velocity ellipsoid is everywhere 
aligned with the confocal ellipsoidal coordinates. When A — > 
—a, the ellipsoidal coordinate surface degenerates into the 
area inside the focal ellipse (Fig.|5Jl. The area outside the fo- 
cal ellipse is labeled by p = —a. Hence, T xx is perpendicular 
to the surface inside and T^^ is perpendicular to the sur- 
face outside the focal ellipse. On the focal ellipse, i.e. when 
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A = fl = —a, both stress components therefore have to be 
equal. Similarly, T MM and T vv are perpendicular to the area 
inside (ju = — jS) and outside (y = —j3) the two branches of 
the focal hyperbola, respectively, and have to be equal on 
the focal hyperbola itself (p = v — —/3). This results in the 
following two continuity conditions 



Txx(—a, -a, v) = T w (-q, -a, u), 



T^{X,-p,-P)=T vv (X,-p,-p). 



(2.17a) 



(2.17b) 



These conditions not only follow from geometrical argu- 
ments, but are also precisely the conditions necessary to 
avoid singularities in the Jeans equations 12.161 when A = 
fi — —a and fj, = v = —p. For the sake of physical under- 
standing, we will also obtain the corresponding continuity 
conditions by geometrical arguments for the limiting cases 
that follow. 



2.5 Limiting cases 

When two or all three of the constants a, j3 or 7 are equal, 
the triaxial Stackel models reduce to limiting cases with 
more symmetry and thus with fewer degrees of freedom. 
We show in H2. 61 that solving the Jeans equations for all the 
models with two degrees of freedom reduces to the same 
two-dimensional problem. EL89 first solved this generalised 
problem and applied it to the disc, oblate and prolate case. 
Evans et al. (2000) showed that the large radii case with 
scale- free DF reduces to the problem solved by EL89. We 
solve the same problem in a different way in and obtain a 
simpler expression than EL89. In order to make application 
of the resulting solution straightforward, and to define a uni- 
fied notation, we first give an overview of the limiting cases. 



2.5.1 Oblate spheroidal coordinates: prolate potentials 

When 7 = f3, the coordinate surfaces for constant A and fi 
reduce to oblate spheroids and hyperboloids of revolution 
around the :r-axis. Since the range of v is zero, it cannot 
be used as a coordinate. The hyperboloids of two sheets 
are now planes containing the a;-axis. We label these planes 
by an azimuthal angle Xi defined as tanx = z/y. In these 
oblate spheroidal coordinates (A, fi, \) the potential Vs has 
the form (cf. Lynden-Bell 1962b) 

/(A) - m 9(X) 



Vs = -- 



(2.18) 



A-/x (A + /3)( M + /3)' 

where the function g(x) is arbitrary, and /(r) = (r+a)G(r), 
with G(t) as in eq. 12.41 . The denominator of the second 
term is proportional to y 2 + z 2 , so that these potentials are 
singular along the entire x-axis unless g(x) = 0. In this case, 
the potential is prolate axisymmetric, and the associated 
density ps is generally prolate as well (de Zeeuw et al. 1986) . 
The Jeans equations 12.161 reduce to 



dTxx Txx 



dX 2(A - fi) 
T\ A 



C^MM _|_ Tp/i 



dfi 



2(,i -A) 



+ 



+ 



Txx - T xx 


dV s 


2(X + P) 




T MM — T xx 


_ dVs_ 


2(p + (3) 


d/2 


dT xx 


dVs 


d\ 





(2.19) 



The continuity condition I2.17al still holds, except that the 
focal ellipse has become a focal circle. For /j, = — /3, the one- 
sheeted hyperboloid degenerates into the i-axis, so that T M/J 
is perpendicular to the a;-axis and coincides with T xx . This 
gives the following two continuity conditions 



Txx(-a,-a,x) = T^{-a, -a,x), 
T^(\,-(3,x) = T xx (\,-p, X )- 



(2.20) 



By integrating along characteristics, Hunter et al. (1990) 
obtained the solution of l|2.19fl for the special prolate models 
in which only the thin I- and O-tube orbits are populated, 
so that T)j M = and Tx\ = 0, respectively (cf. H2.5.61 . 



2.5.2 Prolate spheroidal coordinates: oblate potentials 

When p = a, we cannot use (i as a coordinate and replace 
it by the azimuthal angle <j>, defined as tan<j!> = y/x. Sur- 
faces of constant A and v are confocal prolate spheroids 
and two-sheeted hyperboloids of revolution around the z- 
axis. The prolate spheroidal coordinates (A, (j>, v) follow from 
the oblate spheroidal coordinates (A, /j,, x) by taking fi—>v, 
X — > 4> and P — > a — > 7. The potential Vs(X,^>,f) is (cf. 
Lynden-Bell 1962b) 



Vs = - 



fW - m 



(2.21) 



A — v (X + a)(u + a) 

In this case, the denominator of the second term is propor- 
tional to R — x +y , so that the potential is singular along 
the entire z-axis, unless g(4>) vanishes. When g(4>) = 0, the 
potential is oblate, and the same is generally true for the 
associated density ps- 

The Jeans equations 12.161 reduce to 

dTxx Txx — Tfaf, Txx — T vv _ dVs 
~d\~ 2(A + q) 2(\-v) ~ ^ P ~dX' 



dT* 



dV s 



(2.22) 



dV s 
av 



dTj/y T vv Txx Tw Tfyfy 
dv 2{v — A) 2{v + a) ' 

For A = — a, the prolate spheroidal coordinate surfaces re- 
duce to the part of the z-axis between the foci. The part 
beyond the foci is reached if v = —a. Hence, in this case, 
Txx is perpendicular to part of the z-axis between, and T vv 
is perpendicular to the part of the z-axis beyond the foci. 
They coincide at the foci (A = v — —a), resulting in one 
continuity condition. Two more follow from the fact that 
T^ is perpendicular to the (complete) z-axis, and thus co- 
incides with Txx and T vv on the part between and beyond 
the foci, respectively: 



Txx(—a, 4>, —a) = T v „{— a, 4>, —a) 
T X x{-ct,<t),v) =T^{-a,(p,u), 
T UU (X, 4>, —a) = T^(X,4>,—a). 



(2.23) 



For oblate models with thin S-tube orbits (Txx = 0, see 
H2.5.61 . the analytical solution of 12.221 was derived by 
Bishop (1987) and by de Zeeuw & Hunter (1990). Robijn 
& de Zeeuw (1996) obtained the second-order velocity mo- 
ments for models in which the thin tube orbits were thick- 
ened iteratively. Dejonghe & de Zeeuw (1988, Appendix D) 
found a general solution by integrating along characteristics. 



6 Van de Ven et al. 



Evans (1990) gave an algorithm for solving 112.2211 numeri- 
cally, and Arnold (1995) computed a solution using charac- 
teristics without assuming a separable potential. 

2.5.3 Confocal elliptic coordinates: non-circular discs 

In the principal plane z = 0, the ellipsoidal coordinates re- 
duce to confocal elliptic coordinates (A,/i), with coordinate 
curves that are ellipses (A) and hyperbolae (fi), that share 
their foci on the symmetry y-axis. The potential of the per- 
fect elliptic disc, with its surface density distribution strat- 
ified on concentric ellipses in the plane z — [y — —7), 
is of Stackel form both in and outside this plane. By a su- 
perposition of perfect elliptic discs, one can construct other 
surface densities and corresponding disc potentials that are 
of Stackel form in the plane z = 0, but not necessarily out- 
side it (Evans & de Zeeuw 1992). The expression for the 
potential in the disc is of the form lEHHt with g(x) = 0: 

Vg = _ fW-fM } (224) 

A — Li 

where again f(r) = (r + oi)G(t), so that G(r) equals the 
potential along the y-axis. 

Omitting all terms with v in 12.161 . we obtain the Jeans 
equations for non-circular Stackel discs 



dT X x 


+ 


T\\ — 




dV s 


dX 


2(A- 


-aO 


~ P ^X 


9T^ 


+ 


T MM — 




dVs 


dfx 


2(A»- 


- a) 


" P ^7 



(2.25) 



where now p denotes a surface density. The parts of the y- 
axis between and beyond the foci are labeled by A = —a 
and 11 — —a, resulting in the continuity condition 

T xx (-a, -a) = T MM (-a, -a). (2.26) 




Figure 3. Behaviour of the confocal ellipsoidal coordinates in the 
limit of large radii r. The surfaces of constant X become spheres. 
The hypcrboloids of constant fj, and v approach their asymptotic 
surfaces, and intersect the sphere on the light and dark curves, 
respectively. These form an orthogonal curvilinear coordinate sys- 
tem (fi, v) on the sphere. The black dots indicate the transition 
points (fi = v = — j3 ) between both sets of curves. 



The general Jeans equations in conical coordinates, as de- 
rived by Evans et al. (2000), reduce to 12.291 for vanishing 
mixed second moments. At the transition points between the 
curves of constant fj, and v (fj, = v = —f3), the tensor com- 
ponents TJjp and T vv coincide, resulting in the continuity 
condition 



Txx(r,-0,-P) =7W(r, -0,-/3). 



(2.30) 



2.5.4 Conical coordinates: scale-free triaxial limit 

At large radii, the confocal ellipsoidal coordinates (A, fi, v 
reduce to conical coordinates (r,fi,v), with r the usual dis- 
tance to the origin, i.e., r 2 = x 2 +y 2 + z 2 and jx and v angulai 
coordinates on the sphere (Fig.|2Jl. The potential Vs(r, fi, v 
is scale-free, and of the form 

F{n)-F{u) 



V s = -F(r) + ■ 



■ a)(r - 



(2.27) 

7)G(t), as 



where F(r) is arbitrary, and F(t) = (r 
in eq. (12. 411 . 

The Jeans equations in conical coordinates follow from 
the general triaxial case 12.1611 by going to large radii. Taking 
A — *• r 2 2> —a > fJ-i v i the stress components approach each 
other and we have 



Txx-T vv 1 d__ ]_8_ 

r~* ' dX 2r dX ' 



2(A- M ) ' 2(X-u) (2 ' 28) 

Hence, after multiplying 12.16aH by 2r, the Jeans equations 
for scale-free Stackel models are 





+ 2 rr 






dr 


dr 




r 






dTw + 
dp 


T — 
20* - 


T 

-f) 


dV s 




dT vv 


T-w 


MA* 


dV s 

~ P ^7 




dv 


2(v - 





(2.29) 



2.5.5 One- dimensional limits 

There are several additional limiting cases with more sym- 
metry for which the form of Vs (Lynden-Bell 1962b) and 
the associated Jeans equations follow in a straightforward 
way from the expressions that were given above. We only 
mention spheres and circular discs. 

When a=/3=7, the variables 11 and v loose their mean- 
ing and the ellipsoidal coordinates reduce to spherical coor- 
dinates (r,8,(j)). A steady-state spherical model without a 
preferred axis is invariant under a rotation over the angles 9 
and (f>, so that we are left with only one Jeans equation in r, 
and Tee = T^. This equation can readily be obtained from 
the CBE in spherical coordinates (e.g., Binney & Tremaine 
1987). It also follows as a limit from the Jeans equations 
12.16t for triaxial Stackel models or from any of the above 
two-dimensional limiting cases. Consider for example the 
Jeans equations in conical coordinates 12.291 . and take 
H — > 8 and v — > (f>. The stress components T rr and T MM = 
T vv — T^tf, = Tee depend only r, so that we are left with 

dT rr | 2{T rr -Tee) _ p dV s ^ ^ 

dr r dr ' 

which is the well-known result for non-rotating spherical 
systems (Binney & Tremaine 1987). 

In a similar way, the one Jeans equation for the circular 



General solution Jeans equations 7 



disc-case follows from, e.g., the first equation of (12.2511 by 
taking /j, = —a and replacing T w by T^, where <j> is the 
azimuthal angle defined in H2.5.2I With X+a = R 2 this gives 



cITrb, Trr — _ dVs 
dR R ~ ~ P ~dR' 



(2.32) 



which may be compared with Binney & Tremaine (1987), 
their eq. (4.29). 



2.5.6 Thin tube orbits 

Each of the three tube orbit families in a triaxial Stackel 
model consists of a rotation in one of the ellipsoidal coordi- 
nates and oscillations in the other two ( 4'2.'2t . The I-tubes, 
for example, rotate in v and oscillate in A and /i, with turning 
points fii, fi2 and Ao, so that a typical orbit fills the volume 



"7 < v < —f3, fj.i < (i < fj.2, —at < A < Aq. 



(2.33) 



When we restrict ourselves to infinitesimally thin I-tubes, 
i.e., fii — [i2, there is no motion in the ^i-coordinate. There- 
fore, the second-order velocity moment in this coordinate 
is zero, and thus also the corresponding stress component 
TvL = 0. As a result, eq. 12.16bl reduces to an algebraic 
relation between T\ x and Tl v . This relation can be used to 
eliminate T'„ and T\ x from the remaining Jeans equations 
12.1bal and 12.1(jcl respectively. 

HZ92 solved the resulting two first-order PDEs (their 
Appendix B) and showed that the same result is obtained 
by direct evaluation of the second-order velocity moments, 
using the thin I-tube DF. They derived similar solutions for 
thin O- and S-tubes, for which there is no motion in the 
A-coordinate, so that T xx = and T xx = 0, respectively. 

In Stackel discs we have - besides the flat box orbits 
- only one family of (flat) tube orbits. For infinitesimally 
thin tube orbits Txx = 0, so that the Jeans equations 12. 251 
reduce to two different relations between T MM and the density 
and potential. In >I3.4.4I we show how this places restrictions 
on the form of the density and we give the solution for T M/J . 
We also show that the general solution of 12.2511 . which we 
obtain in Sj3] contains the thin tube result. The same is true 
for the triaxial case: the general solution of 12,161 , which we 
derive in |3J contains the three thin tube orbit solutions as 
special cases C i|4.6.61 . 



2.6 All two-dimensional cases are similar 

EL89 showed that the Jeans equations in oblate and pro- 
late spheroidal coordinates, 12.1911 and 12.2211 . can be trans- 
formed to a system that is equivalent to the two Jeans equa- 
tions 12.251 in confocal elliptic coordinates. Evans et al. 
(2000) arrived at the same two-dimensional form for Stackel 
models with a scale-free DF. We introduce a transformation 
which differs slightly from that of EL89, but has the advan- 
tage that it removes the singular denominators in the Jeans 
equations. 

The Jeans equations 12.191 for prolate potentials can be 
simplified by introducing as dependent variables 



so that the first two equations in 12.191 transform to 
=-(A+/?)*(Ai+/3)* 



dTxx T\\—Tp,p_ 
dX 2(X~fi) 



dT Ul 



T uu —T\, 



dfi 2(fj,-X) 



-(ti+P)Hx+f3y 



dV 3 dT xx 

P dX dX 

dVs | dT xx 

djj, <9/i 



(2.35) 



The third Jeans equation 12.191 can be integrated in a 
straightforward fashion to give the x-dependence of T xx . It 
is trivially satisfied for prolate models with g(x) = 0. Hence 
if, following EL89, we regard T xx (A,/x) as a function which 
can be prescribed, then equations 12.351 have known right 
hand sides, and are therefore of the same form as those of 
the disc case 12.251 . The singular denominator (/x + 0) of 
12.191 has disappeared, and there is a boundary condition 



T w (A,-/3)=0, 



(2.36) 



due to the second continuity condition of 12.201 and the 
definition 12.341 . 

A similar reduction applies for oblate potentials. The 
middle equation of l|2.22|l can be integrated to give the in- 
dependence of T^, and is trivially satisfied for oblate mod- 
els. The remaining two equations 12.221 transform to 



ox 

dT vu 



- + ■ 



2(X-v) 
Tin,— Txx 



=-(A + a) 2 (-a-v) 2 
=-(-a-z/) 5 (A+a)i 



dv ' 2(v-X) 
in terms of the dependent variables 

Tt (A, v) = (X+a)i (-a-v)i (Ttt—a^), 

We now have two boundary conditions 

Txx(-a,v)=0, % v (X,-a)=0, 



dV s 

-dx + 




dX 


dVs 




dv 


dv 


, T = 


- A, v. 



(2.37) 



(2.39) 



as a consequence of the last two continuity conditions of 
<2"35l and the definitions JOfl . 

In the case of a scale-free DF, the stress components in 
the Jeans equations in conical coordinates 12.291 have the 
form T TT — r~''T rT (fi 1 i'), with £ > and r = r,/i,is. After 
substitution and multiplication by r (,+1 , the first equation 
of 12.291 reduces to 



(2 ^T rr -\- Tyfj, -f- T V i/ 



dV s 
or 



(2.40) 



T T t(X, fi) = (X+(3)i(p+P)nT T .~T xx ), 



X,fi, (2.34) 



When C = 2, T rr drops out, so that the relation between T^ 
and T vv is known and the remaining two Jeans equations can 
be readily solved (Evans et al. 2000). In all other cases, T rr 
can be obtained from 12.401 once we have solved the last two 
equations of 12.291 for 7^ M and T vv . This pair of equations 
is identical to the system of Jeans equations 12.251 for the 
case of disc potentials. The latter is the simplest form of the 
equivalent two-dimensional problem for all Stackel models 
with two degrees of freedom. We solve it in the next section. 

Once we have derived the solution of 12.251 . we may 
obtain the solution for prolate Stackel potentials by replac- 
ing all terms —pdVs/dr (r = A, /i) by the right-hand side of 
12.351 and substituting the transformations 12.341 for T\\ 
and T MM . Similarly, our unified notation makes the applica- 
tion of the solution of 12.251 to the oblate case and to models 
with a scale-free DF straightforward ( H3.41 . 
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3 THE TWO-DIMENSIONAL CASE 

We first apply Riemann's method to solve the Jeans equa- 
tions |2.25p in confocal elliptic coordinates for Stackel discs 
( H2.5.31 . This involves finding a Riemann-Green function 
that describes the solution for a source point of stress. The 
full solution is then obtained in compact form by represent- 
ing the known right-hand side terms as a sum of sources. 
In >I3.2I we introduce an alternative approach, the singular 
solution method. Unlike Riemann's method, this can be ex- 
tended to the three-dimensional case, as we show in J2I We 
analyse the choice of the boundary conditions in detail in 
H3.3I In H3.4I we apply the two-dimensional solution to the 
axisymmetric and scale- free limits, and we also consider a 
Stackel disc built with thin tube orbits. 



The combination QCT — TC*Q is a divergence for any twice 
different iable function Q because 



QCT - T£*g = dL/dX + dM/dfi, 



where 



L(X,h) 



M(X,h) 



g or Tdg_ cigr 

2 d/j, 2 dfi A— [i ' 

gdT rag c 2 gr 



(3.7) 



(3.8) 



2 OX 2 dX X-fj, ' 

We now apply the PDE (13.411 and the definition 13.61 in zero- 
subscripted variables Ao and Ho . We integrate the divergence 
13.71 1 over the domain D = {(Ao, Ho)'- A < Ao < oo, H < Ho < 
— a}, with closed boundary P (Fig.^J. It follows by Green's 
theorem that 



3.1 Riemann's method 

After differentiating the first Jeans equation of 12.251 with 
respect to h and eliminating terms in T M/J by applying the 
second equation, we obtain a second-order partial differen- 
tial equation (PDE) for T\\ of the form 



d 2 T X x 



+ 



9T\ 



dXdn 2(X-n) dX 2(A-/x) dfj, 
Here Uxx is a known function given by 
1 d 



Uxx(X,h)- (3.1) 



Uxx = 



(A-At)l dn 



dV s 



2(A-/Lt) dfi 



(3.2) 



We obtain a similar second-order PDE for T MM by inter- 
changing A <-> fi. Both PDEs can be solved by Riemann's 
method. To solve them simultaneously, we define the linear 
second-order differential operator 



C 



if 



Cl 



d 



+ 



(■■2 



d 



dXd/j, X—^dX A— /i <9/i 



(3.3) 



with ci and C2 constants to be specified. Hence, the more 
general second-order PDE 



£T = U, 



(3.4) 



with T and U functions of A and fi alone, reduces to those 
for the two stress components by taking 



T = Txx 



Cl 



3 
2 ' 



r, =§, U = Uxx, 
T — T MM : ci = 4, C2 = §, U = U^. 



(3.5) 



In what follows, we introduce a Riemann-Green function 
g and incorporate the left-hand side of 13.41 into a diver- 
gence. Green's theorem then allows us to rewrite the surface 
integral as a line integral over its closed boundary, which 
can be evaluated if g is chosen suitably. We determine the 
Riemann-Green function g which satisfies the required con- 
ditions, and then construct the solution. 



3.1.1 Application of Riemann's method 

We form a divergence by defining a linear operator £.* , called 
the adjoint of C (e.g., Copson 1975), as 



C* = 



2 



+ 



d f Ci \ d_( c 2 



dXd/i dX\X—nJ dfi\X—fi 



(3.6) 



dA d Mo gC T - TCoG 



dfj, L(X , fio) - <t dXoM(Xo,Ho), (3.9) 



where T is circumnavigated counter-clockwise. Here Co and 
jCq denote the operators 13.31 and 13.61 in zero-subscripted 
variables. We shall seek a Riemann-Green function <5(Ao, Ho) 
which solves the PDE 



CoG = 0, 



(3.10) 



in the interior of D. Then the left-hand side of 13.91 be- 
comes JT dAod/io<5(Ao, Mo) U(Xo, Ho)- The right-hand side 
of 13.91 has a contribution from each of the four sides of the 
rectangular boundary T. We suppose that M(Ao,A*o) and 
L(Ao,/io) decay sufficiently rapidly as Ao — > oo so that the 
contribution from the boundary at Ao = oo vanishes and the 
infinite integration over Ao converges. Partial integration of 
the remaining terms then gives for the boundary integral 



dA 



[(<9Ao Ao— A*o) ] ^° [( 9fM> Xo—fM) ) ] 

A A»o=M fj, A =A 

oo 

X M0 = -Q 

We now impose on g the additional conditions 

Q{X,n) = l, (3.12) 

and 

og c 2 g 



dXo Ao — no 

eg Cl g 



= o 







Ho = H, 



A = A. 



(3.13) 



Oho Xo—Ho 
Then eq. £0D gives the explicit solution 



oo —a 



T(X,h) = JdXo Jdno g(Xo, Ho) U(X ,ho) 



IdXo 

LVc'Ao Ao— Ho 

X MO = -a 



for the stress component, once we have found the Riemann- 
Green function g. 
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Figure 4. The (Xq, fio)-plane. The total stress at a held point 
(A,/i), consists of the weighted contributions from source points 
at (Ao,/io) in the domain D, with boundary T. 



3.1.2 The Riemann-Green function 

Our prescription for the Riemann-Green function G(Xo,p.o) 
is that it satisfies the PDE ffi.lUp as a function of Ao and 
jUo, and that it satisfies the boundary conditions 13.1211 and 
13.131 1 at the specific values Ao = A and /io = /i. Conse- 
quently Q depends on two sets of coordinates. Henceforth, 
we denote it as Q{\, /t; Ao, fio). 

An explicit expression for the Riemann-Green function 
which solves 13.101 is (Copson 1975) 



(Ao-mo) C2 (A- Mo ) c1 " 



(3.15) 



(3.16) 



n i\ \ \ \^u — [iu y/\ — p,v, , 

Q{\H\ A ,Ato) = 7T r- F(w), 

where the parameter w is defined as 

^ _ (Ao-A)(/Jo-p) 
(Ao— )Uo)(A— n) ' 

and F(w) is to be determined. Since w = when Ao = A 
or fj, = n, it follows from 13.121 that the function F has 
to satisfy F(0) = 1. It is straightforward to verify that Q 
satisfies the conditions 13.131 . and that eq. 13.101 reduces to 
the following ordinary differential equation for F(w) 

w(l-w)F" + [l-(2 + ci-c 2 )w]F' - ci{l-c 2 )F = 0. (3.17) 

This is a hypergeometric equation (e.g., Abramowitz & Ste- 
gun 1965), and its unique solution satisfying F(0) = 1 is 



F(w) = 2 F 1 (c 1 ,l-c 2 ;l; 



(3.18) 



The Riemann-Green function l)3.15|l represents the influence 
at a field point at (A,/i) due to a source point at (Ao,Ato). 
Hence it satisfies the PDE 

CQ(X,fi; A ,po) = <5(Ao— \)5(fM>— /j). (3.19) 

The first right-hand side term of the solution 13.141 is a sum 
over the sources in D which are due to the inhomogeneous 
term U in the PDE fTH . That PDE is hyperbolic with char- 
acteristic variables A and ^t. By choosing to apply Green's 
theorem to the domain D, we made it the domain of depen- 
dence (Strauss 1992) of the field point (A,/x) for 113.411 . and 
hence we implicitly decided to integrate that PDE in the 
direction of decreasing A and decreasing fx. 

The second right-hand side term of the solution 13.141 
represents the solution to the homogeneous PDE CT — 
due to the boundary values of T on the part of the bound- 
ary u\ = —a which lies within the domain of dependence. 



There is only one boundary term because we implicitly re- 
quire that T(A,/i) — > as A — > <x. We verify in t|3.1.4l that 
this requirement is indeed satisfied. 

3.1.3 The disc solution 

We obtain the Riemann-Green functions for T\\ and T^,, 
labeled as Q\\ and <5 MM , respectively, from expressions 13.151 
and 13.181 by substitution of the values for the constants 
ci and C2 from l|3.5|l . The hypergeometric function in Q\\ 
is the complete elliptic integral of the second kind 2 , E(w). 
The hypergeometric function in S^m can also be expressed in 
terms of E(w) using eq. (15.2.15) of Abramowitz & Stegun 
(1965), so that we can write 



<5aa(A, fi; Ao, no) 



Q W (A, /x; Ao,/xo) 



(Ap-po)^ 2E(w) 
(A-/t)5 Tr(Ao-At)' 

(Ao-^o)i 2E(w) 
(A-/i)2 7r(A-/i )' 



(3.20a) 



(3.20b) 



Substituting these into 13.141 gives the solution of the stress 
components throughout the disc as 



Tx\(\,fJ-) = 



7r(A— llY- 



°j ^\\ ~H)\pHo 



,s dV s 
aAo 



(Aq-^q)2 dVs ] 

2 p dfi j 



dA 



E(w) 



H = -q 



(A + Q)^[(Ao + a)2T A A(Ao,-a)]|, 

(3.21a) 



oo —a 



7r(A — /X) - 



E(w) 




(A-jUo) 





,s dVs 
-{\o-Ho) 2 p-^ — 

OfM). 



+ 



-ydAo 

X 



E(w) 



(A-Mo) 

fi =—a 



^-[(A +a)*T w (A 0: 



(A -/t )2 dVs \ 

2 p aA j 

-e,)\}. (3.21b) 



This solution depends on p and Vs, which are assumed 
to be known, and on T\\(\,— a) and T MM (A, — a), i.e., the 
stress components on the part of the y-axis beyond the 
foci. Because these two stress components satisfy the first 
Jeans equation of 12.251 at /i = —a, we are only free to 
choose one of them, say T MM (A, — a). T\\(\, — a) then fol- 
lows by integrating this first Jeans equation with respect to 
A, using the continuity condition 12.261 and requiring that 
T\\(X, —a) — » as A — > oo. 



3.1.4 Consistency check 

We now investigate the behaviour of our solutions at large 
distances and verify that our working hypothesis concerning 
the radial fall-off of the functions L and M in eq. 13.81 is 



We use the definition E{w) = / 2 d6 y 1 — w sin 2 
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correct. The solution 13.141 consists of two components: an 
area integral due to the inhomogeneous right-hand side term 
of the PDE (I3.4|> . and a single integral due to the boundary 
values. We examine them in turn to obtain the conditions 
for the integrals to converge. Next, we parameterise the be- 
haviour of the density and potential at large distances and 
apply it to the solution 13.211 and to the energy equation 
12.1U1 to check if the convergence conditions are satisfied for 
physical potential-density pairs. 

As Ao — > oo, w tends to the finite limit (po — p)/(X — p). 
Hence E(w) is finite, and so, by (13.2UI . Qxx = 0(X^ 2 ) and 
<5 M n = C(Ao /2 ). Suppose now that Uxx(X ,po) = 0(X^ h ~ 1 ) 
and ?7 MM (Ao, po) = 0(X^ mi ~ 1 ) as Ao — > oo. The area in- 
tegrals in the solution 13.1411 then converge, provided that 
h > 7j and mi > |. These requirements place restrictions 
on the behaviour of the density p and potential Vs which 
we examine below. Since Gxx(X, p; Ao, po) is 0(A~ 1,/2 ) as 
A — > oo, the area integral component of T\a(A, p) behaves as 



0{\- 1/2 f™ Ao !l " 1/2 dA ) and so is 0(A~ !l ). Similarly, with 
S(j/i(A, p; Ao, po) = 0(A~ 3 / 2 ) as A — > oo, the first component 
oiT^{X,p) isO(A -" 11 ). 

To analyse the second component of the solution 13.141 . 
we suppose that the boundary value 7aa(Ao, —a) = 0(Aq~' 2 ) 
and r,j M (Ao, —a) = C(Aq m2 ) as Ao — + oo. A similar analysis 
then shows that the boundary integrals converge, provided 
that h > \ and 7712 > § , and that the second components of 
T X x(X,p) and T MM (A, p) are 0{X~ h ) and 0(X~ m2 ) as A -> 
00, respectively. 

We conclude that the convergence of the integrals in the 
solution H3.140 requires that Txx(X,p) and T MM (A, p) decay 
at large distance as 0(X~ l ) with I > ± and 0(A~ m ) with 
m > I, respectively. The requirements which we have im- 
posed on U (Ao, (Mo) and T(Ao, —a) cause the contributions to 
<f r dpoL(Xo, po) in Green's formula (13.911 from the segment 
of the path at large Ao to be negligible in all cases. 

Having obtained the requirements for the Riemann- 
Green function analysis to be valid, we now investigate the 
circumstances in which they apply. Following Arnold et al. 
(1994), we consider densities p that decay as N(p)X~ s ^ 2 at 
large distances. We suppose that the function G(r) intro- 
duced in eq. 12.41 is 0(t 6 ) for — i < 5 < as r — > 00. 
The lower limit 5 — — i corresponds to a potential due to a 
finite total mass, while the upper limit restricts it to poten- 
tials that decay to zero at large distances. 

For the disc potential 12.241 . we then have that f(r) = 
0(t s+1 ) when r — > 00. Using the definition 13.211 . we obtain 



Uxx(X,p) 



f'M-f'W + Vs + f'W dp 



2{\-pY 



(X-fi) dp' 



/'(A) -/'(/*). Vs + f'(p) dp 



2{X-pY 



(A dX' 



(3.22a) 



(3.22b) 



where p is the surface density of the disc. It follows that 
Ux\(X,p) is generally the larger and is 0(X S ~ S ^ 2 ~ 1 ) as 
A —* 00, whereas f/ Mfl (A,/i) is 0(A _2_S / 2 ). Hence, for the 
components of the stresses 13.211 we have Txx = 0(X S ~ S ^ 2 ) 
and Tfj,^ = C(A -1 ~ S / 2 ). This estimate for Uxx assumes that 
dp/dp is also 0(X~ S ^ 2 ). It is too high if the density becomes 
independent of angle at large distances, as it does for discs 
with s < 3 (Evans & de Zeeuw 1992) . Using these estimates 
with the requirements for integral convergence that were 



obtained earlier, we obtain the conditions s > 25 + 1 and 
s > 1, respectively, for inhomogeneous terms in Txx(X, p) 
and T M/J (A, p) to be valid solutions. The second condition 
implies the first because S < 0. 

With V s (X,p) = 0(X S ) at large A, it follows from the 
energy equation 12,101 for bound orbits that the second- 
order velocity moments (v 2 ) cannot exceed 0(X S ), and hence 
that stresses T TT = p{v 2 ) cannot exceed 0(X S ~ S ^ 2 ). This im- 
plies for Txx(X,p) that s > 25 + 1, and for r,j M (A,/i) we have 
the more stringent requirement that s > 28 + 3. This last 
requirement is unnecessarily restrictive, but an alternative 
form of the solution is needed to do better. Since that al- 
ternative form arises naturally with the singular solution 
method, we return to this issue in H3.2.6I 

Thus, for the Riemann-Green solution to apply, we find 
the conditions s > 1 and — | < S < 0. These conditions are 
satisfied for the perfect elliptic disk (s = 3, 8 = — |), and for 
many other separable discs (Evans & de Zeeuw 1992). 



3.1.5 Relation to the EL89 analysis 

EL89 solve for the difference A = Txx — Tfifi using a Green's 
function method which is essentially equivalent to the ap- 
proach used here. EL89 give the Fourier transform of their 
Green's function, but do not invert it. We give the Riemann- 
Green function for A in Appendix A, and then rederive it by 
a Laplace transform analysis. Our Laplace transform anal- 
ysis can be recast in terms of Fourier transforms. When we 
do this, we obtain a result which differs from that of EL89. 



3.2 Singular Solution Superposition 

We have solved the disc problem 12.251 by combining the two 
Jeans equations into a single second-order PDE in one of the 
stress components, and then applying Riemann's method to 
it. However, Riemann's method and other standard tech- 
niques do not carry over to a single third-order PDE in one 
dependent variable, which is the best that one could expect 
to have in the general case. We therefore introduce an alter- 
native but equivalent method of solution, also based on the 
superposition of source points. In constrast to Riemann's 
method, this singular solution method is applicable to the 
general case of triaxial Stackel models. 



3.2.1 Simplified Jeans equations 
We define new independent variables 
Sxx{X,p) = \X — p\ ' Taa(A, p), 
S^{X,p) = \p-X\z T^(X,p), 



(3.23) 



where |.| denotes absolute value, introduced to make the 
square root single- valued with respect to cyclic permutation 
of A — > p — > A. The Jeans equations 12.251 can then be 
written in the form 



as. 



Su 



x li dV s N 
dX 2(X-p) =-|A-H 2 P^ > -=3i(A,M), 



dS u 



Sx 



dp 2(p-X) 



(3.24a) 



(3.24b) 
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For given density and potential, gi and <?2 are known func- 
tions of A and p. Next, we consider a simplified form of 
l|3.24fl by taking for gi and 172, respectively 



§i(X,li) = 0, ff 2 (A,/tt) = 8(\ -X)5(no-n), 



(3.25) 



with — p < (i < (Mo < —a < A < Ao. A similar set of simpli- 
fied equations is obtained by interchanging the expressions 
for gi and <?2- We refer to solutions of these simplified Jeans 
equations as singular solutions. 

Singular solutions can be interpreted as contributions 
to the stresses at a fixed point (A, pi) due to a source point 
in (Ao,/xo) (Fig. [IJ. The full stress at the field point can be 
obtained by adding all source point contributions, each with 
a weight that depends on the local density and potential. In 
what follows, we derive the singular solutions, and then use 
this superposition principle to construct the solution for the 
Stackel discs in >I3.2.6I 

3.2.2 Homogeneous boundary problem 

The choice 113.2511 places constraints on the functional form 
of S\\ and S^m- The presence of the delta-functions in (72 
requires that S^p contains a term — <5(Ao— X)H.(ij,o — n), with 
the step- function 



H(x — xo) 



0, x < xo, 

1, x > Xo- 



(3.26) 



Since Tt'(y) = 5(y), it follows that, by taking the partial 
derivative of — <5(Ao — X)TL(iio — n) with respect to [i, the 
delta-functions are balanced. There is no balance when S\\ 
contains <5(Ao — A), and similarly neither stress components 
can contain <5(/io — t 1 )- We can, however, add a function of A 
and n to both components, multiplied by H(\q—\)H(hq—i£). 
In this way, we obtain a singular solution of the form 



(3.27) 



in terms of functions A and B that have to be determined. 
Substituting these forms in the simplified Jeans equations 
and matching terms yields two homogeneous equations 

dA B „ dB A 
9A ~ 2{X-n) 



0, 



and two boundary conditions 
1 



0, 



(3.28) 



A(\o,v) 



B(A,mo) =0. 



(3.29) 



2(A -/Lt)' 

Two alternative boundary conditions which are useful below 
can be found as follows. Integrating the first of the equations 
l|3.28|l with respect to A on 11 = fio, where B(X, no) = 0, gives 



A{X,no) 



1 



(3.30) 



2(A -^o) ' 

Similarly, integrating the second of equations l|3,28|l with 
respect to n on A = Ao where A is known gives 

Mo - A* 



S(Ao,/x) = 



(3.31) 



4(A -a*o)(Ao-m) ' 
Even though expressions I3.3UH and 13. 3111 do not add new 
information, they will be useful for identifying contour inte- 
gral formulas in the analysis which follows. 



plane 




Figure 5. Contours C M and C A in the complex z-plane which 
appear in the solution 13.371) . The two cuts running from p, to no 
and one from A to Aq make the integrands single-valued. 



We have reduced the problem of solving the Jeans equa- 
tions l)2.25|l for Stackel discs to a two-dimensional bound- 
ary problem. We solve this problem by first deriving a one- 
parameter particular solution ( >I3.2.31 and then making a 
linear combination of particular solutions with different val- 
ues of their free parameter, such that the four boundary 
expressions are satisfied simultaneously ( H3.2.4I . This gives 
the solution of the homogeneous boundary problem. 

3.2.3 Particular solution 

To find a particular solution of the homogeneous equations 
ll3~25t with one free parameter z, we take as an Ansatz 



A(X,n) cc (X-n) ai (z-\r(z-n) a3 , 

B(X,n) oc (\-fi) bl (z~V b2 (z-ri b3 , 

with <2i and bi (i = 1, 2, 3) all constants. Hence, 



(3.32) 



dA _ f ax a 2 

dX \X—fx z—X 

dB = /_6i &3_ 

dn \fj,— X z—fj, 



1 



2(A-m) 
1 



2aiA 



2bxB 



Z — fJ, 

z-X 
z-X 



(3.33) 



2(/i— A) \ z— /i 

where we have set 02 = — Oi and 63 = — &!• Taking ai 
61 = |, the homogeneous equations are satisfied if 



z-X A 



(z-xy 



-b 2 



(3.34) 



so, a 3 



B p (X,n) = 



•|. We denote the resulting solutions as 



(2 — A) 2 (z—fi) 5 
\fi— X\i 



(3.35a) 



(3.35b) 



(z—fj,) 2 (z—X) 2 

These particular solutions follow from each other by cyclic 
permutation A — > n — > A, as is required from the symmetry 
of the homogeneous equations l|3.28^ . 

3.2.4 The homogeneous solution 

We now consider a linear combination of the particular solu- 
tion 13.35H by integrating it over the free parameter z, which 
we assume to be complex. We choose the integration con- 
tours in the complex z-plane, such that the four boundary 
expressions can be satisfied simultaneously. 
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We multiply B p (X,ii) by (z — jUo) 5 , and integrate it 
over the closed contour (Fig.|5Jl. When fi = hq, the inte- 
grand is analytic within C M , so that the integral vanishes by 
Cauchy's theorem. Since both the multiplication factor and 
the integration are independent of A and /t, it follows from 
the superposition principle that the homogeneous equations 
are still satisfied. In this way, the second of the boundary 
expressions 13.291 is satisfied. 

Next, we also multiply B P (X, fx) by (z — Ao) - 5 , so that 
the contour C A (Fig. |5J encloses a double pole when A = 
Ao- From the Residue theorem (e.g., Conway 1973), it then 
follows that 



-B (Xo,n)dz = 

(z-XoP 



27Tl(Ao — /i) ; 



_d_ 

dl 



Z — fXQ 

z — (1 



(z — fl )2 (X — fl)2 

(z-fi)^(z-Xo) 2 

m(fxo-fx) 
(A -/to) 3 (Ao-/i) 



dz = 



(3.36) 



which equals the boundary expression 13.311 . up to the fac- 
tor 47ri(Ao— /io) 2 • 

Taking into account the latter factor, and the ratio 
13.341 of A and B, we postulate as homogeneous solution 



A(X,») = 



4iri 



(z-fi ) 2 dz 



I A — Aio 1 2 (z-X)2(z-n)*(z-X y- 



B(X,fi) = 



l/i- A|i 



4ni \x - 



/'() 



(2-/10) 2 dz 



(«-/t)3(«-A)3(*-Ao) ? 



(3.37a) 



(3.37b) 



with the choice for the contour C still to be specified. 

The integrands in 13.371 consist of multi-valued func- 
tions that all come in pairs (z— t) 1 ^ 2 "" 1 (z~ To) 1 / 2- ™, for in- 
teger m and n, and for r being either A or fi. Hence, we can 
make the integrands single- valued by specifying two cuts in 
the complex z- plane, one from fi to fio and one from A to Ao. 
The integrands are now analytic in the cut plane away from 
its cuts and behave as z~ 2 at large distances, so that the 
integral over a circular contour with infinite radius is zero 3 . 
Connecting the simple contours C A and with this cir- 
cular contour shows that the cumulative contribution from 
each of these contours cancels. As a consequence, every time 
we integrate over the contour C A , we will obtain the same 
result by integrating over — C M instead. This means we inte- 
grate over C* and take the negative of the result or, equally, 
integrate over in clockwise direction. 

For example, we obtained the boundary expression for 
B in 13.361 by applying the Residue theorem to the double 
pole enclosed by the contour C A . The evaluation of the inte- 
gral becomes less straightforward when we consider the con- 
tour — C M instead. Wrapping the contour around the branch 
points fi and /to (Fig.|3Jl, one may easily verify that the con- 
tribution from the two arcs vanishes if their radius goes to 
zero. Taking into account the change in phase when going 
around the two branch points, one may show that the con- 
tributions from the two remaining parts of the contour, par- 
allel to the real axis, are equivalent. Hence, we arrive at the 



3 We evaluate the square roots as (z— t) ■ 
with | arg(z — t)| < 7r. 



\z— t\ exp i arg(z — r) 



Im 



c 



z -plane 



Re 



T 



T 



/ 



Figure 6. Integration along the contour C T . The contour is 
wrapped around the branch points r and to (t = A, u.), and split 
into four parts. Ti and T3 run parallel to the real axis in opposite 
directions. T2 and T4 are two arcs around r and to, respectively. 



following (real) integral 



S(A ,/i) 



1 (A-Mo) 2 



dt 



2 ^ (Ao-Mo)* J i^o-t) 2 V t-n 



Ho—t 



The substitution 



t = n + 



Guo-/t)(Ao-/to) sin 2 



(3.38) 



(3.39) 



(no— /u) s— (Ao— m) 

then indeed gives the correct boundary expression l)3.31|l . 

When we take fi — fio in I3.37bl . we are left with the 
integrand (z — A)~ 3//2 (z — Ao)~ 1//2 . This is analytic within 
the contour C M and hence it follows from Cauchy's theorem 
that there is no contribution. However, if we take the con- 
tour — C x instead, it is not clear at once that the integral 
indeed is zero. To evaluate the complex integral we wrap 
the contour C A around the branch points A and Ao (Fig.|SJ. 
There will be no contribution from the arc around Ao if its 
radius goes to zero. However, since the integrand involves 
the term z — A with power — |, the contribution from the 
arc around A is of the order e~ 1//2 and hence goes to in- 
finity if its radius e > reduces to zero. If we let the two 
remaining straight parts of the contour run from A + e to Ao, 
then their cumulative contribution becomes proportional to 
tanS(e), with 9(e) approaching ^ when e reduces to zero. 
Hence, both the latter contribution and the contribution 
from the arc around A approaches infinity. However, care- 
ful investigation of their limiting behaviour shows that they 
cancel when e reaches zero, as is required for the boundary 
expression f?(A,/io) = 0. 

We have shown that the use of C A and — C M gives the 
same result, but the effort to evaluate the contour integral 
varies between the two choices. The boundary expressions 
for A(X,fi), 13.291 and 13.301 are obtained most easily if we 
consider C A when A = Ao and — C M when fi — [io- In both 
cases the integrand in 13.37al has a single pole within the 
chosen contour, so that the boundary expressions follow by 
straightforward application of the Residue theorem. 

We now have proven that the homogeneous solution 
13.371 solves the homogeneous equations 13.281 . satisfies the 
boundary values 13.291 - 13.311 separately and, from the ob- 
servation that C A and — C M produce the same result, also 
simultaneously. 
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3.2.5 Evaluation of the homogeneous solution 

The homogeneous solution 13.371 consists of complex con- 
tour integrals, which we transform to real integrals by wrap- 
ping the contours C x and C M around the corresponding pair 
of branch points (Fig. |SJ . To have no contribution from the 
arcs around the branch points, we choose the (combination 
of) contours such that the terms in the integrand involv- 
ing these branch points have powers larger than —1. In this 
way, we can always evaluate the complex integral as a (real) 
integral running from one branch point to the other. 

In the homogeneous solution I3.37al for A we choose 
C = C x and in Ij3.37b|l for B we take C = -C"\ Taking 
into account the changes in phase when going around the 
branch points, we obtain the following expressions for the 
homogeneous solution 



A(A,m) = 



1 |A-/i|* 



dt 



t-no 



2T|Ao-Mo|5i*-MV(*->0(*-M)(*o-t) : 



B(A,/*) = 



1 MO 

A-^|2 f dt 



Ho—t 



2 * |Ao-uojWA-tV(A-t)(t-M)(Ao-i)' 



(3.40a) 



(3.40b) 



By a parameterisation of the form 13.391 . or by using an 
integral table (e.g., Byrd & Friedman 1971), expressions 
13.401 can be written conveniently in terms of the complete 
elliptic integral of the second kind, E, and its derivative E' 



A(X,n;X ,no) 



B(X,/j,; X ,fio) = 



E(w) 

7r(Ao— fJ.) ' 

2wE'(w) 
'tt(Ao-A)' 



(3.41a) 



(3.41b) 



with w defined as in 13.161 . The second set of arguments 
that were added to A and B make explicit the position 
(Xo,Ho) of the source point which is causing the stresses at 
the field point (A,/x). 



3.2.6 The disc solution 

The solution of equations 13.241 with right hand sides of the 
simplified form 



<?i(A,m) = <5(A -A)5Gu -m)> ffa(A,|i) = °> 



(3.42) 



is obtained from the solution 13.271 by interchanging A <-» fi 
and Ao <-> [lo- It is 

S\\=B(n, A; fjbo, Xo)H(Xo—X)H(fJ.o~Li) — 8(lJ-o—lJ-)'H(Xo—X), 
S^Ain, A; no, X )H(Xo-X)H(ii -n). ( 3 ' 43 ) 

To find the solution to the full equations 13.241 at (A, /j,), 
we multiply the singular solutions 13. 2711 and 13.4311 by 
<7i(Ao,/io) and #2(Ao,/io) respectively and integrate over D, 
the domain of dependence of (A,/j). This gives the first two 
lines of the two equations 13.441 below. The terms in the 
third lines are due to the boundary values of S 1 ,^ at fi = —a. 
They are found by multiplying the singular solution 13.271 
evaluated for fio — —a by — S^^Ao, — a) and integrating 
over Ao in D. It is easily verified that this procedure correctly 
represents the boundary values with singular solutions. The 



final result for the general solution of the Jeans equations 
13.241 1 for Stackel discs, after using the evaluations 13.411 1. is 



S\x(X,fj,) = — dX gi(Ao,M) 



+ JdX Jdfj.o 

X M 



.2wE'(w) , E(w) 

-gi(Xo, Mo) — 7 r+ff2(A , 



dAo S^(Xo, —a) 



tt(mo-m) 
E(w) 



7r(A -/x)_ 

H0=-ot 



w(Xo-n) 
(3.44a) 



S W (A, (i) — - 

oo —a 

+ JdXoJdno 

A j-i 



— OL 

j d/xo 



ff2(A,/Uo) 



,. E(w) ,2wE'(w) 
-5i(Ao,Mo)— 7T r-32(A ,Mo' 



+ S A i M (A, —a) — ydAo S MA ,(Ao, —a) 



7f(A — ho) ' 7r(Ao — A) 

2wE'(w) 



tt(Ao-A) 

fl =-a 



(3.44b) 



The terms (/io~ m) -1 an d (Ao— A) -1 do not cause singularities 
because they are canceled by components of w. In order to 
show that equations 13.4411 are equivalent to the solution 
13.211 1 given by Riemann's method, integrate the terms in 
E'(w) by parts, and use the definitions of S TT , gi and g2- 

3.2.7 Convergence of the disc solution 

We now return to the convergence issues first discussed 
in H3.1.4I where we assumed that the density p decays as 
iV(/i)A~ S//2 at large distances and the Stackel potential as 
0(X 5 ). For the physical reasons given there, the assigned 
boundary stress T MM (A, — a) cannot exceed 0(X S ~ S ^ 2 ) at 
large A, giving an S^(X,~a) of 0{X s ' s/2+1/2 ). It follows 
that the infinite integrals in S l /1M (Ao, — a) in the solution 
13.4411 require only that s > 28 + 1 for their convergence. 
This is the less restrictive result to which we referred earlier. 

The terms in the boundary stress are seen to contribute 
terms of the correct order 0(X s ~ s/2+1/2 ) to Sxx{X,fi) and 
S7i/j(A,p). The formulas for the density and potential show 
that gi(X,fl) = 0{X s - s/2 - 1/2 ) while g 2 (A,/Li) is larger and 
£)^-s/2-i/2-j ag ^ _^ oo. The Ao integrations with g\ and g2 
in their integrands all converge provided s > 28 + 1. Hence, 
both Saa(A, h) and S Mfl (A, fj.) are 0{X s - 3/2+1/2 ), so that the 
stress components T TT (X, fx) (r = A, (j) are 0(X S ~ S ^ 2 ), which 
is consistent with the physical reasoning of H3.1.4I 

Hence, all the conditions necessary for 13.441 to be a 
valid solution of the Jeans equations 13.241 for a Stackel 
disc are satisfied provided that s > 28 + 1. We have seen in 
H3.1.4l that 8 must lie in the range [— \, 0). When 5^0 the 
models approach the isothermal disk, for which also s — 1 
when the density is consistent with the potential. Only then 
our requirement s > 28 + 1 is violated. 



3.3 Alternative boundary conditions 

We now derive the alternative form of the general disc so- 
lution when the boundary conditions are not specified on 
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= — ol but on /i = —f3, or on A = —a rather than in the 
limit A — > oo. While the former switch is straightforward, 
the latter is non-trivial, and leads to non-physical solutions. 



3. 3. 1 Boundary condition for fi, 

The analysis in H3.1l and >!3.2l is that needed when the bound- 
ary conditions are imposed at large A and at [i = —a. The 
Jeans equations 12. 2511 can be solved in a similar way when 
one or both of those conditions are imposed instead at the 
opposite boundaries A = —a and/or /i = — f3. The solution 
by Riemann's method is accomplished by applying Green's 
theorem to a different domain, for example D' = {(Ao,/Uo): 
A < Ao < oo, — f3 < fio < fJ,} when the boundary conditions 
are at (X = — j3 and as A -— > oo. The Riemann-Green func- 
tions have to satisfy the same PDE iTHrJI and the same 
boundary conditions 13.121 and 13.131 . and so again are 
given by equations 13.2Ual and 13.20bl . The variable w is 
negative in D' instead of positive as in D, but this is unim- 
portant. The only significant difference in the solution of 
eq. 13.41 is that of a sign due to changes in the limits of the 
line integrals. The final result, in place of eq. 13.141 . is 



T(X,fj.) = — JdXoJdfioG(Xo,(M))U(Xo,fJ-o) 



-0 oo 



9T c 2 T 

i • 

A M>=— P 



+ 



A{) — (AO 



(3.45) 



To apply the method of singular solutions to solve for the 
stresses when the boundary stresses are specified at fj, = — j3 
rather than at /i = —a, we modify the singular solutions 
13.271 by replacing the step- function Tt(fio— A*) by —H(^—fj,o) 
throughout. No other change is needed because both func- 
tions give — 5(n — jUo) on partial differentiation with respect 
to fi. The two-dimensional problem for A and B remains 
the same, and so, as with Riemann's method, its solution 
remains the same. Summing over sources in D 1 now gives 



S\\{X,n) = ~y d ^o 3i(A ,m) 

A 

oo n 

dX d/i Q 



.2wE'(w) , E(w) 

5i(A ,Moj— 7 r +P2(A ,mo 



A -p 

OO 

ydA S MM (Ao, — fi) 



tt(^o-m) 
E(w) 



■k{X —h) 

M)=— 



7r(Ao-M)_ 

(3.46a) 



S w (A,p) = Jdfi g 2 (X,Ho) 

-0 



oo n 

JdXojdfjLo 



.. E{w) .2wE'{w) 
-ffi(Ao,/io)-7T r - 52(A ,/jo) 



+ S Mi (\,-/3)- /dA 5 MM (A 0j -/3) 



as an alternative to equations 13.441 . 



7r(A — (Ho) 7f(Ao— A) 

2wE'(w) 



tt(Ao-A) 

M0 = 



(3.46b) 



3.3.2 Boundary condition for X 

There is a much more significant difference when one assigns 
boundary values at A = —a rather than at A — * oo. It is still 
necessary that stresses decay to zero at large distances. The 
stresses induced by arbitrary boundary data at the finite 
boundary A = —a do decay to zero as a consequence of 
geometric divergence. The issue is that of the rate of this 
decay. We find that it is generally less than that required by 
our analysis in H3.1.4I 

To isolate the effect of boundary data at A = —a, 
we study solutions of the two-dimensional Jeans equations 
12.251 when the inhomogeneous right hand side terms are set 
to zero and homogeneous boundary conditions of zero stress 
are applied at either fj, = — a or (i = - (3. These solutions 
can be derived either by Riemann's method or by singular 
solutions. The solution of the homogeneous PDE CT = is 



T(A,,)=-/d, [(g-^)e(A,,;Ao,Mo)] , (3.47) 



A = -a 



for the case of zero stress at = —a, and 



T (W^(f^^) e(w °'H ' (3 - 48) 



for the case of zero stress at n = — f3. 

The behaviour of the stresses at large distances is gov- 
erned by the behaviour of the Riemann-Green functions Q 
for distant field points (A, fi) and source points at Ao = —a. 
It follows from equations ||33H that T xx {X,fi) = 0(A~ 1/2 ) 
and T Mfl (A,^i) = 0(A~ 3//2 ). As a restult, the radial stresses 
dominate at large distances and they decay as only the in- 
verse first power of distance. Their rate of decay is less than 
0(X S ~ S / 2 ) - obtained in H3.1.4l from physical arguments - if 
the requirement s > 25+ 1 is satisfied. This inequality is the 
necessary condition which we derived in >!3.2.6l for fElH to 
be a valid solution of the disc Jeans equations 13.241 . It is 
violated in the isothermal limit. 

There is a physical implication of radial stresses which 
decay as only the inverse first power of distance. It implies 
that net forces of finite magnitude are needed at an outer 
boundary to maintain the system, the finite magnitudes aris- 
ing from the product of the decaying radial stresses and the 
increasing length of the boundary over which they act. That 
length grows as the first power of distance. Because this sit- 
uation is perhaps more naturally understood in three dimen- 
sions, we return to it in our discussion of oblate models in 
H3.4.2I For now, lacking any physical reason for allowing a 
stellar system to have such an external constraint, we con- 
clude that boundary conditions can be applied only at large 
A and not at A = —a. 



3.3.3 Disc solution for a general finite region 

We now apply the singular solution method to solve 
equations 13.241 in some rectangle ^t m i n < H < fJ>m,ax> 
A m in < A < A max , when the stress 5 MM is given a boundary 
in /i, and Saa is given on a boundary in A. This solution 
includes 13.441 and 13.461 as special cases. It will be needed 
for the large-radii scale-free case of H3.4.3I 
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As we saw in >I3.3.1I singular solutions can easily be 
adapted to alternative choices for the domain of dependence 
of a field point (X,h)- Originally this was D, the first of the 
four quadrants into which (Ao, /Lto)-space is split by the lines 
Ao = A and Ho = h (Fig. It has the singular solution 
13.271 . We then obtained the singular solution for the fourth 
quadrant D' simply by replacing H(ho — h) by —7i{fi — fio) 
in (13.2711 . We can similarly find the singular solution for 
the second quadrant A m i n < Ao < A, fj, < Ho < /imax by 
replacing H.(Xo—X) by — H(X— Ao), and for the third quadrant 
A m in < A < A, ^ min < Ho < H by replacing H(Xo~X) by 
— 7i(X— Ao) and7i.(fio—l-i) by —Ti,(fi—fio)- We find the part of 
the solution of equations 13.241 due to the right hand side 
g terms by multiplying the first and second terms of the 
singular solutions by gi(Xo,Ho) and g 2 (Xo, Ho), respectively, 
and integrating over the relevant domain. We use A = A e 
and /J, — He to denote the boundaries at which stresses are 
specified. We find the part of the solution generated by the 
boundary values of S M(i by multiplying the singular solution 
13.271 1. modified for the domain and evaluated at ho — fi e , 
by ±5 , M(J (Ao, He) and integrating over Ao in the domain. The 
plus sign is needed when /i e = /i m i n and the minus when 
He = Mmax- Similarly, the part of the solution generated 
by the boundary values of S\\ is obtained by multiplying 
the singular solution 13.431 . modified for the domain and 
evaluated at Ao = A e , by ±S\\(X e , ho) and integrating over 
Ho in the domain. The sign is plus if A e = A m i n and minus 
if A e = A max . The final solution is 



^e 

S\\(\,fi) = S\\(\ e ,Li) — JdXogi(Xo, fj-) 

A 

A e fj, e 

+ Jd\o Jdfio [gi(Xo,fj,o)B(fj,,X;fj.o,Xo)+g2(Xo,fJ,o)A(X,fi;Xo,fj,o)] 



3.4 Applying the disc solution to limiting cases 



A f-i 
Ae 



dXoS^ fl (Xo,^ e )A(X,H;Xo,fJ,e)~JdHoSx\(X e ,fj,o)B(H,X;fio,X e ), 

" (3.49a) 

S^(X,n) = Sw(X,iJ, e ) - Jdii g 2 (X, no) 

A e Me 

+ JdX Q jdfj.o [gi(Xo,Ho)A(H,X;Ho,Xo)+g2(Xo,Ho)B(X,H;Xo,Ho)] 

A j-i 

Ae Me 

—JdX S fll _ l (Xo,He)B(X,H;Xo,He)—JdHoS\\(X e ,Ho) A(H,X;Ho,X r )- 
x M (3.49b) 

This solution is uniquely determined once gi and g2 are 
given, and the boundary values S M(i (Ao, He) and S\\{X e , no) 
are prescribed. It shows that the hyperbolic equations 13.241 
can equally well be integrated in either direction in the 
characteristic variables A and h- Solutions 13.441 and 13.461 
are obtained by taking A e — > oo, S\\(X e , Ho) — * 0, setting 
He. = ~ot and /i e = — /3 respectively, and evaluating A and 
B by equations 13.411 . 



We showed in §221 that the Jeans equations for prolate and 
oblate potentials and for three-dimensional Stackel models 
with a scale-free DF all reduce to a set of two equations 
equivalent to those for the Stackel disc. Here we apply our so- 
lution for the Stackel disc to these special three-dimensional 
cases, with particular attention to the behaviour at large 
radii and the boundary conditions. This provides further in- 
sight in some of the previously published solutions. We also 
consider the case of a Stackel disc built with thin tube orbits. 

3.4-1 Prolate potentials 

We can apply the disc solution 13.461 to solve the Jeans 
equations 12.351 by setting S\\{X,h) — |A — h\ 5 ^AA(A,/i) 
and S fl ^(X, ft) = \h — A| 27^ M (A, h), and taking 



gi(X,H)=-\X-^(X+(3)2(H+P)- 



g 2 (X, h)=-\h~M 2 (X+(3)-2(h+P)- 



dVs dT x 



+ - 



dX 
dVs_ 

&h ' <9/i 



+ 



dX 



(3.50) 



The boundary terms in 5^ (A, — (5) vanish because of the 
boundary condition 12.361 . As before, we regard the az- 
imuthal stress T xx as a variable that can be arbitrarily as- 
signed, provided that it has the correct behaviour at large A 
f H3.1.41 . The choice of T xx is also restricted by the require- 
ment that the resulting solutions for the stresses T\\ and 
T MM must be non-negative (see H2.31 . 

The analysis needed to show that the solution obtained 
in this way is valid requires only minor modifications of 
that of ^3.2.71 We suppose that the prescribed azimuthal 
stresses also decay as 0{X s ~ a ^ 2 ) as A — » oo. As a result 
of the extra factor in the definitions 13.5U1 . we now have 
gi(X,n) = 0(X s - s/2 ) and <? 2 (A, M ) = 0{X' 3 ' 2 ) as A -> oo. 
The Ao integrations converge provided s > 28 + 2, and 5*aa 



and S u 



s/2+l 



). Hence the stresses T\\ and T M 



are 0(X 

which follow from T TT = T xx + S TT / yj '(A- /i)(A + P)(h+P), 
are once again 0{X s ~ a / 2 ). The requirement s > 25 + 2 is 
no stronger than the requirement s > 28 + 1 of H3.2.7I it is 
simply the three-dimensional version of that requirement. It 
also does not break down until the isothermal limit. That 
limit is still 5^0, but now s — *• 2. 



3.4-2 Oblate potentials 

The oblate case with Jeans equations 12.371 differs sig- 
nificantly from the prolate case. Now S\\{X, v) = A — 
H 5 ^aa(A, i^) vanishes at A = —a and S^A,^) = \u — 
X\^T VV (X, v) vanishes at v = —a. If one again supposes that 
the azimuthal stresses can be assigned initially, then 
one encounters the problem discussed in >!3.3.2l of excessively 
large radial stresses at large distances. To relate that anal- 
ysis to the present case, we use the solution 13.441 with h 
replaced by v, and with zero boundary value 5^^(A, — a), 
and for gi and g2 the right hand side of 12.371 multiplied by 
A — v\ 2 and \v — A| 5 , respectively. 

The estimates we obtained for the prolate case are 
still valid, so the stresses T AA and T vv are C(A 4 " s/2 ). Dif- 
ficulties arise when this solution for S\ \ does not van- 
ish at A = —a, but instead has some nonzero value k(v) 



16 Van de Ven et at 



there. To obtain a physically acceptable solution, we must 
add to it a solution of the homogeneous equations 12.371 
with boundary values T\\(— a,v) — —k(v)/s/—(X — v and 
T vv (A, —a) — 0. This is precisely the problem we discussed 
in H3.3.2l where we showed that the resulting solution gives 
T\x (A, p) = 0(\- 1/2 ), and hence T xx (\, fi) = C(A _1 ). This 
is larger than 0(X S ~ 3 ^ 2 ) when the three-dimensional require- 
ment s > 28 + 2 is met. We therefore conclude that the ap- 
proach in which one first selects the azimuthal stress 
and then calculates the other two stresses will be unsuc- 
cessful unless the choice of is fortunate, and leads to 
k(u) = 0. Otherwise, it leads only to models which either 
violate the continuity condition T\\ — = at A = —a, 
or else have radial stresses which require external forces at 
large distances. 

The physical implication of radial stresses which decay 
as only C^A" 1 ), or the inverse second power of distance, is 
that net forces of finite magnitude are needed at an outer 
boundary to maintain the system. This finite magnitude 
arises from the product of the decaying radial stresses and 
the increasing surface area of the boundary over which they 
act, which grows as the second power of distance. This sit- 
uation is analogous to that of an isothermal sphere, as il- 
lustrated in problem 4-9 of Binney & Tremaine (1987), for 
which the contribution from an outer surface integral must 
be taken into account in the balance between energies re- 
quired by the virial theorem. 

There are, of course, many physical models which sat- 
isfy the continuity condition and whose radial stresses decay 
in the physically correct manner at large distances, but some 
strategy other than that of assigning initially is needed 
to find them. In fact, only Evans (1990) used the approach 
of assigning initially. He computed a numerical solution 
for a mass model with s = 3 and Vs oc 0(A -1 / 2 In A) for 
large A, so that the stresses there should be C(A _2 lnA). 
He set = — ipVs, which is of this magnitude, and inte- 
grated from A = —a in the direction of increasing A for a 
finite range. Evans does not report on the large A behaviour, 
and it is possible that his choice of gives k(v) = 0, but 
his Figure 2 especially shows velocity ellipsoids which be- 
come increasingly elongated in the radial direction, consis- 
tent with our prediction that T\\ generally grows as C(A _1 ) 
when the boundary value of T\\ is assigned at A = — a. 

A more common and effective approach to solve the 
Jeans equations for oblate models has been to specify the 
ratio T\\/Tm, and then to solve for one of those stresses 
and Tfaj, (Bacon, Simien & Monnet 1983; Dejonghe & de 
Zeeuw 1988; Evans & Lynden-Bell 1991; Arnold 1995). This 
leads to a much simpler mathematical problem with just 
a single first-order PDE. The characteristics of that PDE 
have non-negative slopes dX/dv, and therefore cut across the 
coordinate lines of constant A and v. The solution is obtained 
by integrating in along the characteristics from large A. The 
continuity conditions 12,231 are taken care of automatically, 
the region —7 < v < — a < 00 is covered, and it is easy to 
verify that the stresses so obtained are everywhere positive. 



3.4-3 Large radii limit with scale-free DF 

We found in H2.5.4I that the first of the Jeans equations in 
conical coordinates (12.2911 reduces to an algebraic relation 
for the radial stress T rr . The problem that remains is that 



of solving the second and third Jeans equations for T M(J and 
T vv . Those equations are exactly the same as those of the 
disc case after we apply the coordinate permutation A — > 
p, — > v, and the physical domain is —7 < v < —fi < p < 
—a with finite ranges of both variables. Hence, the solution 
13.491 can be applied with T MM assigned at either p e = —a 
or p e = —fi, and T vv at either v £ — —fi or v a — —7. For 
<7i and gi we take the same expressions as for the disc case, 
i.e., the right-hand side of 13.241 . but with A — > p — > v and 
multiplied by r^ . To obtain T M(J and T„„ from the S\\ and 
respectively, we use the transformation 



(3.51) 



with C, > the scaling factor. We can choose to specify 
the stress components on the two boundaries /1 = —fi and 
v — —fi. For a given radius r these boundaries cover the 
circular cross section with the (x, z)-plane (Fig.^J. We can 
consider the (x, z)-plane as the starting space for the solu- 
tion. It turns out that the latter also applies to the triax- 
ial solution ( H4.6.31 and compares well with Schwarzschild 
(1993), who used the same plane to start his numerically 
calculated orbits from. 



3.4-4 Thin tube orbits 

For infinitesimally thin tube orbits in Stackel discs we have 
that S\\ = f H2.5.6H . so that equations 13.241 reduce to 



-^-p ffl (A, M ), ^ =fla (A, M ). 



(3.52) 



A solution is possible only if the right hand side terms satisfy 
the subsidiary equation 

9a (A, M ) = -2^[(A-^i(A,/i)]. (3.53) 

We find below that this equation places restrictions on the 
form of the (surface) density p, and we use this relation 
between gi and 52 to show that the disc solution 13.441 yields 
the right results for the stress components. 

If we write the disc potential 12.241 as a divided differ- 
ence, Vs = —/[A, fi], we have that 

gi = (A-/i)5p/[A,A,/Lt], g 2 = (\- (i)? pf[\, p]. (3.54) 

Upon substitution of these expressions in 13.531 we obtain a 
PDE in p, of which the solution implies the following form 
for the density 



p(X,fl) 



/(A) 



(3.55) 



where /(A) is an arbitrary function independent of p. 
From 13.521 and the definition 13.231 it then follows that 



T^(X,p, v) 



-2/(A)v7[A, X,p\. The tube density that 



de Zeeuw, Hunter & Schwarzschild (1987) derive from the 
DF for thin tube orbits in the perfect elliptic disk (their 
eq. [4.25]) is indeed of the form 13.551 . 

To show that the general disc solution 13.441 gives 
S\\(X, p) — 0, we substitute eq. I|3.53|l for 32 (A, p) in l|3.44a|l . 
After partial integration and using 



2(Ao — po) 



d E{w) _ 2wE'(w) 



dpo 7r(A — p) ^{po — p)' 



(3.56) 
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we find that the area integral reduces to 



oo 

JdXoi.gi(X , ^) - 2(A +a)ffi(A 0) -a) 



E(w) 



tt(Xo—h) 

fj, =—a 



(3.57) 



The first part cancels the first line of 13.44al and since from 
13.521 we have that — 2(Ao + a)gi(Xa, — a) = S MM (Ao, — a), 
the second part cancels the third line. Hence, we have 
Sxx(X,fi) = as required. To see that the general disc so- 
lution also yields Sf l f 1 (X,n) correctly, we apply similar steps 
to 13.44bll . where we use the relation 

d 2wE'(w) 



-2(A — Ho): 



E(w) 



duo 7r(A — A) 7r(A— jito) 
We are finally left with 



(3.58) 



— at 

S' MM (A, /i) = 5 MM (A, -a) - J d/i O g 2 (A,^ ), (3.59) 

which is just the second equation of 13.521 integrated with 
respect to fi. 



4 THE GENERAL CASE 

We now solve the system of three Jeans equations 12.161 
for triaxial Stackel models by applying the singular solu- 
tion superposition method, introduced in A3. 21 for the two- 
dimensional case. Although the calculations are more com- 
plex for a triaxial model, the step-wise solution method is 
similar to that in two dimensions. Specifically, we first sim- 
plify the Jeans equations and show that they reduce to 
a three-dimensional homogeneous boundary problem. We 
then find a two-parameter particular solution and apply con- 
tour integration to both complex parameters to obtain the 
general homogeneous solution. The latter yields the three 
singular solutions of the simplified Jeans equations, from 
which, by superposition, we construct the general solution. 

4.1 Simplified Jeans equations 

We start by introducing the functions 



SW(A,/i, v) = y/ (X-h)(X-u)(h-u) T TT (A, fx, v), 



with r = A, /i, v, to write the Jeans equations for triaxial 
Stackel models l|2.16[l in the more convenient form 



dSxx 




Si/ is 


ox 


2(A-/x) 


2{X-v) 




Suu 


Sxx 


dfi 


2(n~v) 


2(A*-A) 


dS vv 


S\\ 


Sfj.fi 


dv 


2{v-X) 


2{v-n) 



= gi(x,v, v), 



= g2(x,n,u), 



= g3(X,n,u), 



where the function <?i is defined as 



gi(X,,u,v) = - V (A-/x)(A-i/)(/i-i/) p 



dVs 

ox ' 



(4.2a) 



(4.2b) 



(4.2c) 



(4.3) 



and Q2 and gz follow by cyclic permutation A — > /i — » v — > A. 
We keep the three terms A — /x, X—v and fj,— v under one square 
root. With each cyclic permutation two of the three terms 



change sign, so that the combination of the three terms is 
always positive real. Therefore, the square root of the com- 
bination is always single- valued, whereas in the case of three 
separate square roots we would have a multi- valued function. 

We simplify equations 114.211 by substituting for gi, g^ 
and gz, respectively 

gi(X,H,v) = 0, 

g 2 {X,n,v) = S(X -X)S(p. -fJ.)S(v ~v), (4.4) 
gz(X,ti,u) = 0, 

with 

-7 < v < vq < -/3 < h < fio < -a < X < X . (4.5) 

We obtain two similar systems of simplified equations by 
cyclic permutation of the left-hand side of 14.21 . Once we 
have obtained the singular solutions of the simplified system 
with the right-hand side given by 14.41 . those for the other 
two systems follow via cyclic permutation. 

4.2 Homogeneous boundary problem 

The choice 14.41 implies that the functions S TT (X, fi, v) 14.11 
must have the following forms 

Sxx = A{X, 1 i,v)H(X -X)H(v,q-ii)H(vq-v) 

+ F(X,fj,)5(v -u)n(Xo-X)'H(Ho- l u.), 
= B(X,h,v)H{X -X)H(ho-h)'H(vo-v) 

+ G(X,fM)8(v -p)H(Xo-X)n(fi - l u) 

+ H(fi,u)5(Xo-X)H(fio-n)H(vo-u) 

- 6(X - X)8(v — v)Tl{no- v), 
S vv = C(X, fi,v)H{Xo — X)H(fio— n)H(v — v) 

+ I(p.,v)5{X -X)H{p.o-p)TL{v -v), 

with A, B, C and F, G, H, I yet unknown functions of 
three and two coordinates, respectively, and TL the step- 
function 13.261 . After substituting these forms into the sim- 
plified Jeans equations and matching terms we obtain 14 
equations. Eight of them comprise the following two homo- 
geneous systems with two boundary conditions each 



(4.6) 



dF 



(4- 1 ) 7» 



dG 



and 



dH 



01 



G 


2(A- 


/<) 


F 




20*- 


X) 


I 




2(M- 


v) 


H 





= 0, 



2(Ao-M)' 



(4.7) 



= 0, 



0, 



G(A,mo) = 0, 

H(p , v) = 0, 
v ) 



1 



(4.8) 



, dv 2(i> — fu) 2{va — fu) 

We have shown in how to solve these two-dimensional 
homogeneous boundary problems in terms of the complete 
elliptic integral of the second kind E and its derivative E' . 
The solutions arc 



F(X, fi) = 
H{n,v) = 



E(w) 



7r(A() - fl) ' 

2uE'(u) 
tt(vo — v) ' 



G(A, M ) = - 



I(fi,v) 



2wE'(w) 
7r(Ao — A) ' 

E(u) 



(4.9) 
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where u and similarly v, which we will encounter later on, 
follow from w 13.161 1 by cyclic permutation A — > n — > v — > A 
and Ao —* flo vo —* Ao, so that 



(/jo — v){vq — v) 



{vg-v)(\ -\) 



(4.10) 



0*o— fo)0* — ") ' (Ao— Md)(A— i/)' 

The remaining six equations form a three-dimensional homo- 
geneous boundary problem, consisting of three homogeneous 
Jeans equations 



dA 



dB 



dC 



B 




C 




2(A- 


/<•) 


2(A- 


v) 


C 




A 




2(p- 


") 


20*- 


A) 






B 




2(j/- 


A) 


2{y- 


/<) 



= o, 
= o, 

= 0. 



(4.11) 



and three boundary conditions, specifically the values of 
A(Xo, fi,v), B(X, iM),u), and C(A,/i, fo). As in H3.2.2I it is 
useful to supplement these boundary conditions with the val- 
ues of A, B, and C at the other boundary surfaces. These are 
obtained by integrating the pairs of equations 14.1111 which 
apply at those surfaces, and using the boundary conditions. 
This results in the following nine boundary values 



A(X ,fi, v) 
A{X,iio,v) 
A(X,fi,vo) 
B(X ,fJ., v) 
B(X,(j, , v) 
B(X,n,v ) 

C(Xa,n,v) 
C(X,fi , v) 
C(X,n,vo) 



1 

2^ 

J_ 

2^ 



E{u) 



- + - 



2uE'(u) 



(Xo-v)(n-v ) (Xo-h)(vq-v) 
E(v) 2vE'(v) 



{X q -v)(iiq-vq) + (X -fMi)(vo-v) 

Ao— 



E(w) 


X—/i 


4tt(Ao-m) 


_(A- 


^o)(p- 


fo) 


uE'{u) 










(A - 


-Mo)(Ao 


-m) 


0, 








wE'(w) 




Mo-M 




2tt(A -A) 


(/xo- 


-vo)(ji- 


-fo) 


E(u) 








&ir(fj,—i/o) 


(Ao- 


-ju)(A - 





!/ — V 



(4.12) 



Ao-A 



+ - 



Mo-fo 



1 

2^ 
1 

2^ 



2vE'(v) 



(Ao— /io)(A— fo) (jUo— fo)(Ao— A) 
E(w) 2wE'(w) 



(Ao— m)(A — fo) (/*— ^o)(Ao— A) 



If we can solve the three homogeneous equations 14. lit and 
satisfy the nine boundary expressions l|4.12^ simultaneously, 
then we obtain the singular solutions (14.611 . By superposi- 
tion, we can then construct the solution of the Jeans equa- 
tions for triaxial Stackel models. 



4.3 Particular solution 

By analogy with the two-dimensional case, we look for par- 
ticular solutions of the homogeneous equations 14. lit and by 
superposition of these particular solutions we try to satisfy 
the boundary expressions 14.121 simultaneously, in order to 
obtain the homogeneous solution for A, B and C. 



4-3. 1 One-parameter particular solution 
By substitution one can verify that 



A*^u) = £ ^- V ^- U \ t C«-A) (4 . 13) 
(A-/x)(A-z/) {z-fj,){z-v) 

with B p and C p following from A p by cyclic permutation, 
solves the homogeneous equations 14. Ill . To satisfy the nine 
boundary expressions 14.121 , we could integrate this partic- 
ular solution over its free parameter z, in the complex plane. 
From H3.2.4I it follows that, at the boundaries, this results in 
simple polynomials in (X,n,v) and (Ao, fJ-o, vo). This means 
that the nine boundary expressions 14.121 cannot be satis- 
fied, since in addition to these simple polynomials they also 
contain E and E' . The latter are functions of one variable, 
so that at least one extra freedom is necessary. Hence, we 
look for a particular solution with two free parameters. 



4-3.2 Two-parameter particular solution 

A particular solution with two free parameters Z\ and zi 
can be found by splitting the ^-dependent terms of the one- 
parameter solution 14.131 into two similar parts and then 
relabelling them. The result is the following two-parameter 
particular solution 



a p = \/(A-m)(A-^)(m-^) -q 



(zi-xy- 



(A-mXA-!/) 



b p = V , (A-/*)(A-t/)(/j-^) -q 



i=i (z;-m) 2 (zi-vY 
{zi-fj,)i 



0i-i/)(//-A) 



c p _ \/(A-m)(A-^)(m-^) tt 



i=i (zi-u)2(zi-X)i 

2 



(4.14) 



-!/)2 



{v-X){v-ii) 



1 (z i — X)2(z i ~fl) i - 



These functions are cyclic in (A, fi, u), as is required from the 
symmetry of the homogeneous equations 14.111 . The pres- 
ence of the square roots, such as occurred earlier in the solu- 
tion 13.321 for the disc case, allows us to fit boundary values 
that contain elliptic integrals. 

To show that this particular solution solves the ho- 
mogeneous Jeans equations, we calculate the derivative of 
A p (X,fj,, v) with respect to A: 



dA* 



2 



+ 



1 



dX 2 \ A— zi A — z-2 
This can be written as 



A — fi X — v 



(4.15) 



dA p 
dX 



2(A-/x) 
1 



+ 



2{X-v) 



(z 1 -fl)(z 2 -tl){X-u) A P 

(z 1 -X)(z 2 -X)(ij,-v) 

(z 1 -u)(z 2 -v)(X-n) a p 
(z 1 -X)(z 2 -X)(fi-v) 



(4.16) 



From the two-parameter particular solution we have 



A p 
A p 



(£i - y)(z 2 - n)(X - v) 
[z\ — X)(Z2 — X)(p — v) ' 

(z-i - v)(z 2 — v)(X — fi) 
(Z\ — X){Z2 — X)(fl — v) ' 



(4.17) 



so that, after substitution of these ratios, the first homoge- 
neous equation of 14.111 . is indeed satisfied. The remaining 
two homogeneous equations can be checked in the same way. 
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4.4 The homogeneous solution 

In order to satisfy the four boundary expressions of the two- 
dimensional case, we multiplied the one-parameter partic- 
ular solution by terms depending on Ao, /io and the free 
complex parameter 2, followed by contour integration over 
the latter. Similarly, in the triaxial case we multiply the two- 
parameter particular solution 13. 3511 by terms depending on 
Ao, Ato, and the two free parameters zx and Zi, in such a 
way that by contour integration over the latter two complex 
parameters the nine boundary expressions 14. 121 can be sat- 
isfied. Since these terms and the integration are independent 
of A, /i and v, it follows from the superposition principle that 
the homogeneous equations 14.111 remain satisfied. 

The contour integrations over zi and 22 are mutually 
independent, since we can separate the two-parameter par- 
ticular solution 14.1411 with respect to these two parameters. 
This allows us to choose a pair of contours, one contour 
in the zi-plane and the other contour in the 22-plane, and 
integrate over them separately. We consider the same sim- 
ple contours as in the disk case (Fig. around the pairs 
of branch points (A, Ao) and (fi, po), and a similar contour 
around (y, vq). We denote these contours by Cf, Cf and 
C" respectively, with i = 1, 2 indicating in which of the two 
complex planes we apply the contour integration. 



4-4-1 Boundary expressions for B 

It follows from 14.121 that B — at the boundary /i = fiQ. 
From Cauchy's theorem, B would indeed vanish if, in this 
case, in either the zi-plane or 22-plane the integrand for 
B is analytic within the chosen integration contour. The 
boundary expression for B at v — vq follows from the one 
at A = Ao by taking v <-> A and pq <-> Ao. In addition to 
this symmetry, also the form of both boundary expressions 
puts constraints on the solution for B. The boundary ex- 
pressions can be separated in two parts, one involving the 
complete elliptic integral E' and the other consisting of a 
two-component polynomial in r and tq (t = A,/i,^). Each 
of the two parts follows from a contour integration in one of 
the two complex planes. For either of the complex parame- 
ters, z\ or 22 , the integrands will consist of a combination of 
the six terms Zi — r and Zi — To with powers that are half-odd 
integers, i.e., the integrals are of /ij/perelliptic form. If two 
of the six terms cancel on one of the boundaries, we will be 
left with an elliptic integral. We expect the polynomial to 
result from applying the Residue theorem to a double pole, 
as this would involve a first derivative and hence give two 
components. This leads to the following Ansatz 



B(X,n,u) oc - x 

(fi-u)(fi-X) 

(zi—fi)? (21-A0) 5 dzi 



(ssi — u)' (zi-A)2 (zi-^o) 2 {zx — voY- 



(z 2 — jtt) 2 (Z2 — fo) 2 d22 



C-2 



(Z2-f) 2 (22 -A) 2 (Z2-(M)) 2 (22-A0) 2 



(4.18) 



Upon substitution of /1 = /10, the terms involving fio cancel 
in both integrals, so that the integrands are analytic in both 
contours Cf and C% ■ Hence, by choosing either of these 



contours as integration contour, the boundary expression 
B(\,fio, v) = is satisfied. 

When A = Ao, the terms with Ao in the first inte- 
gral in 14.181 cancel, while in the second integral we have 
(z2 — Ao) -2 . The first integral is analytic within Ci, so that 
there is no contribution from this contour. However, the in- 
tegral over Cf is elliptic and can be evaluated in terms of 
E' (cf. H3.2.5H . We apply the Residue theorem to the second 
integral, for which there is a double pole inside the contour 
C2 • Considering Cf 1 and C2 as a pair of contours, the ex- 
pression for B at A = Ao becomes 



B(\,jj,, v) oc — 16ty' 



\j (Ao — /io) (Ao — vq) (ho —vo) 
(^o-^o)(mo-A ) 



uE'{u) 


Mo— ^ 




(Ao— /io)(Ao— At) 



u — V 



(4.19) 



which is the required boundary expression up to a scaling 
factor. As before, we keep the terms Ao — /*o, Ao — ^0 and 
flo—fo under one square root, so that it is single- valued with 
respect to cyclic permutation in these coordinates. 

The boundary expression for B at v — vq is symmetric 
with the one at A = Ao, so that a similar approach can 
be used. In this case, for the second integral, there is no 
contribution from C2 , whereas it can be expressed in terms 
of E' if C2 = C*2 • The first integrand has a double pole in 
Ci . The total contribution from the pair (C" ,0%) gives the 
correct boundary expression, up to a scaling factor that is 
the same as in (14.1911 . 

Taking into account the latter scaling factor, this shows 
that the Ansatz 14.181 1 for B produces the correct boundary 
expressions and hence we postulate it as the homogeneous 
solution for B. The expressions for A and C then follow 
from the ratios 14.17H . Absorbing the minus sign in 14.191 
into the pair of contours, i.e., either of the two contours we 
integrate in clockwise direction, we postulate the following 
homogeneous solution 



A(X fi v)= it^Zl^it^zM / ( a -m)(A-^)(m~^) x 
167r 2 (A— a«) (A— v) Y (Ao— /io)(A — ^o)(mo— ^0) 

(2i-A)2(2i-Ao)2 dzi 



Ci 



C'2 



(zi— At) 2 (21 -f) 2 (zx-fM>) 2 (zx-vo) 2 

(22 — A) 2 (22— tA)) 2 d22 
(z2-fl)2(Z2-v)i (Z2~ fMo)i (Z2- X )i 



(4.20) 



S / A v \ = (Ato-fo)(AtQ-Ao) / {X-n){X-v){n-v) ^ 
16n 2 (xi — — A) y (Ao— Ato)(Ao— ^o)(Ato— va) 

(2i-At)^(zi-A )5 dzx 



Ci 



c 2 



(zx-v) 2 (21 -A) 2 (zx—fio) 2 {zx—uo) 2 

(22— At) 2 (22 — ^0) 2 d2 2 

(22 -f) 2 (22 -A) 2 (22- Ato) 2 {z2-X )i 



(4.21) 
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16n 2 (v- X)(f— n) y (Ao— /tto)(Ao— ^o)(Mo— ^o) 
(zi-i/)i(zi-A )i dzi 



Ci 



C'2 



— A) 2 (jzi-^i) 2 (zi — /xo) 2 (zi — vo) 2 

(Z2~ v)l (z 2 — Vp)^ dz 2 

11 13 

(Z2-A) 2 ( 22 -/i) 2 (z 2 -Ato)2 (Z2-A())2 



(4.22) 



The integrands consist of multi-valued functions that all 
come in pairs of the form (z — r) 2 ~ m ( z — to) 2 _n , for integers 
m and n, with r equal to A, /i or 1/. Hence, completely analo- 
gous to our procedure in >I3.2,4I we can make the integrands 
single-valued by specifying, in the complex zi-plane and Z2- 
plane, three cuts running between the three pairs (A,Ao), 
(fx, Ho), (y,vo) of branch points, that are enclosed by the 
simple contours. The integrands are now analytic in the cut 
plane away from its cuts and behave again as z~ 2 at large 
distances, so that the integral over a circular contour with 
radius going to infinity, will be zero. Hence, connecting the 
simple contours , C? and C" with this circular contour, 
shows that their cumulative contribution cancels 



0. 



1,2. 



(4.23) 



This relation will allow us to make a combination of con- 
tours, so that the nine boundary expressions 114. 1211 can be 
satisfied simultaneously f H4.4.3H . Before doing so, we first 
establish whether, with the homogeneous solution for A and 
C given by (I4.2UI and 14.221 , respectively, we indeed satisfy 
their corresponding boundary expressions separately. 



combination Cf does not give the correct boundary ex- 
pression. We again split both integrals to obtain the required 
complete elliptic integrals. In the first we substitute 



Ao — vo , s Ao — fio, % 
zi - A = (zi— fi ) (zi — vo). 



Hq — vo 



fio — vo 



(4.26) 



For the contour Ci , the first integral after the split can be 
evaluated in terms of E'(v), but the second integral we leave 
unchanged. For the integral in the z 2 -plane we substitute 



Z'l 



Ao — vq. fio — vo , . \ 
vo = -; — (Z2 — fJ-o) — T (22 — Ao). 



Ao — Ho 



Ao — fJ-o 



(4.27) 



We take C2 as contour, and evaluate the first integral after 
the split in terms of E(v). We again leave the second integral 
unchanged. Except for the contour choice, it is of the same 
form as the integral we left unchanged in the zi-plane. 

To obtain the required boundary expression for A at 
fl = fio, it turns out that we have to add the contribution 
of three pairs of contours, C£(7£, C^Cl and C%. With 
the above substitutions 14.261 and I4.27H . the first two pairs 
together provide the required boundary expression, but in 
addition we have two similar contour integrals 



i/87T 



(2 — A) 2 dz 



( Ao- vo) 2 ( A — v) 2 J (z — v) 2 (z — Ao) 2 (z — vo) 2 (2 - fio) 



(4.28) 



with contours C A and C , respectively. The third pair, 
C1C2 , involves the product of two single pole contributions. 
The resulting polynomial 



4-4-2 Boundary expressions for A and C 

The boundary expressions and the homogeneous solution of 
C, follow from those of A by taking A <-» v and Ao <-> vo- 
Henceforth, once we have checked the boundary expressions 
for A, those for C can be checked in a similar way. 

Upon substitution of A = Ao in the expression for A 
14.201 1. the first integrand is proportional to 21 — A' and thus 
is analytic within the contour Ci . The contribution to the 
boundary expression therefore needs to come from either 
or Ci . The substitution 



2ni (A — fio) - 



\a — v t v Ao — fi , \ 
21 - Ao = (zi — fi) (zi — v), 



(4.24) 



fi — v [l~ v 
splits the first integral into two complete elliptic integrals 



(Ao-^o) 2 (A-j/) 2 (fio-v) 2(A -/io) 2 (no-v ) 2 



(4.29) 



can be written in the same form as 14.28H . with contour C M . 
As a result, we now have the same integral over all three 
contours, so that from I4.23H . the cumulative result vanishes 
and we are left with the required boundary expression. 

The expression for A at v = vo resembles the one for B 
at the same boundary. This is expected since their boundary 
expressions in I4.12H are also very similar. The first integral 
now has a contribution from a double pole in the contour 
Ci*. The second integral has no contribution from the con- 
tour C2. However, within C% , the second integral can be 
evaluated in terms of E(w). We obtain the correct bound- 
ary expression A(X,)i,vo) by considering the pair — CfC^ 1 ■ 



Ao — v 



(21 -fi) 2 dzi 



A* v J (z!-v)2 (zi-/i ) 2 {zi-v y. 

Ci 



(z!-v)z dzi 



Ap-/i 

t x ~ u J (zi-/x)4(zi-/xo)*(zi-!4))* 



(4.25) 



Within the contour , the integrals can be evaluated in 
terms of E'(u) and E(u) respectively. When A = Ao, the 
second integral in 14.2UH has a single pole contribution from 
the contour C2. Together, — CJ'C^ , exactly reproduces the 
boundary expression A(Xo,fi, v) in l|4.12|l . 

When fi = fj-Oi both integrands in the expression for 
A have a single pole within the contour Cf . However, the 



4-4-3 Combination of contours 

In the previous paragraphs we have constructed a homoge- 
neous solution for A, B and C, and we have shown that 
with this solution all nine boundary expressions can be sat- 
isfied. For each boundary expression separately, we have de- 
termined the required pair of contours and also contours 
from which there is no contribution. Now we have to find 
the right combination of all these contours to fit the bound- 
ary expressions simultaneously. 

We first summarise the required and non-contributing 
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pairs of contours per boundary expression 



A(X ,fj,, v) 
A(\,fj, ,v) 
A(X,fj,,u ) 



X O2 =t C-i O2 , 
1^2 31 u l °2 > 



— L»j O2 =C O2 , 
! 2 ± Oi O2 , 



B(Xo,n, v, 

B(X, no, v) : ±Cf C 2 T ± CTC^, (4.30) 
B(\,n, u ) 

C{X ,y, v) 
C(X,fi ,u) 
C(X,n,vo) 

where r can be A, /i or v. At each boundary separately, 
X — Xq, ji — no and ^ = vo, the allowed combination of con- 
tours matches between A, B and C. This leaves the question 
how to relate the combination of contours at the different 
boundaries relative to each other. 

From 14. 231 . we know that in both the complex zi-plane 
and 22-plane, the cumulative contribution of the three sim- 
ple contours cancels. As a consequence, each of the following 
three combinations of integration contours 



CT C* 2 M = - CI ( C2 + C2 ) = - ( ci + cr ) C£ 



(4.31) 



will give the same result. Similarly, we can add to each com- 
bination the pairs C\C 2 and C^C 2 , to obtain 

C^C^+CiC^+C^C^C^C^-C^C^C^Ci-C^C^. (4.32) 

The first combination of contour pairs matches the allowed 
range for /i = /in. m 14.301 and the second and third match 
the boundary expressions A = Ao and v — vq. This completes 
the proof that the expressions 114. 21)1 - 14.221 1 for A, B and 
C solve the homogeneous equations 14.111 and satisfy the 
nine boundary expressions 14.121 simultaneously when the 
integration contour is any of the three combinations 14.321 . 
We shall see below that the first of these combinations is 
preferred in numerical evaluations. 



4.5 Evaluation of the homogeneous solutions 

We write the complex contour integrals in the homogeneous 
solutions A, B and C 14.2UH4.22l as real integrals. The re- 
sulting complete hyperelliptic integrals are expressed as sin- 
gle quadratures, which can be evaluated numerically in a 
straightforward way. We also express the complete elliptic 
integrals in the two-dimensional homogeneous solutions F, 
G, H and I 14.91 in this way to facilitate their numerical 
evaluation. 



4-5.1 From complex to real integrals 

To transform the complex contour integrals in 14.201 - 14.221 
in real integrals we wrap the contours C A , C^ and C around 
the corresponding pair of branch points (Fig. |BJ. The inte- 
grands consist of terms z — r and z — To, all with powers 
larger than —1, except z\ — vo and z 2 — Ao, both of which 
occur to the power — ~. This means that for all simple con- 
tours CI (t = X,fj,,v;i = 1,2), except for C\ and C 2 , the 
contribution from the arcs around the branch points van- 
ishes. In the latter case, we are left with the parts parallel 



to the real axis, so that we can rewrite the complex inte- 
grals as real integrals with the branch points as integration 
limits. The only combination of contours of the three given 
in 14.321 that does not involve both C\ and C 2 , is 



S = C? Ca + C?C£ + Cf C% 



(4.33) 



We have to be careful with the changes in phase when wrap- 
ping each of the simple contours around the branch points. 
One can verify that the phase changes per contour are the 
same for all three the homogeneous solutions A, B and C, 
and also that the contribution from the parts parallel to the 
real axis is equivalent. This gives a factor 2 per contour and 
thus a factor 4 for the combination of contour pairs in S. 
In this way, we can transform the double complex contour 
integration into the following combination of real integrals 

-V) MO MO ' y MO MO 

J I d-zidza = ±{j dti Jdt 2 + J dti J dt ^-J d*i J dfa), (4.34) 

5 A M M v M M 

with ti the real part of Zi. 

We apply this transformation to 14.201 - 14.221 . and we 
absorb the factor of 4 left in the denominators into the in- 
tegrals, so that we can write 

A(X,ii,v;Xq , M o ,vo) = { ^° ){ ^' Xo) !' (A 1 A 2 +A s A 4 - A 2 A 3 ) , 
da \ -\ (mo— ^o)(mo— Ao)A 

B(A,^,i^;A ,Mo,t'o)=^7 — ^-(B- L B 2 +B 3 B4-B 2 B 3 ), 

C(X,h,u;X , M o ,u ) = {fM r " o) jg= A "\ A (d CMC, - C 2 C 3 ) , 



■K 2 {y — X)(y — y) 



(4.35) 



where Ai, Bi and d (i = 1, 2, 3, 4) are complete hyperelliptic 
integrals, for which we give expressions below, and 



A 



(x-n)(x-v)(n-v) 

( Ao— mo) (Ao— ^o) (^o— va) 



(4.36) 



The second set of arguments added to A, B and C make 
explicit the position (Ao, fJ-o, vq) of the source point which is 
causing the stresses at the field point (A,/i, v). 



4-5.2 The complete hyperelliptic integrals 

With the transformation described in the previous section 
the expression for, e.g., the complete hyperelliptic integral 
A 2 is of the form 



A, 



MO 

i r At 



(X-t)(t-v ) 



2JX ~t \ (n -t)(t-fj,)(Xo-t)(t-v) 



(4.37) 



The integrand has two singularities, one at the lower inte- 
gration limit t — ji and one at the upper integration limit 
t = fM)- The substitution t = fi + (/io — fi) cos 2 8 removes 

:2(a»o - At)d0. 
Bi and d (i = 



both singularities, since cXtj^J (fiQ—t)(t — fj.) 

All complete hyperelliptic integrals Ai 
1,2,3,4) in H4.35fl are of the form l|4.37|l and have at most 
two singularities at either of the integration limits. Hence, we 
can apply a similar substitution to remove the singularities. 
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This results in the following expressions 



Ai = (X -X) 



tt/2 

2 f sin 2 9 cos 2 9d8 
13A1 



tt/2 



tt/2 



Aa=(vo-v) 



z 2 sin 2 6A6 

2lA z 



A 2 = 



yly4.de 
2/3 Ay 



tt/2 
i 2 



(4.38a) 



yi^y 



tt/2 



S 1 = (A -A) 



B4 = (v — V 



X2 sin 



B 2 = {h 



tt/2 



2/i cos 2 #d6l 
2/3 A H 



t/2 



24 sin OdO 

2lA z 



tt/2 

S3 = (/X0-M)|^ 



(4.38b) 



2/3 cos 

2/1 A y 



tt/2 



d = (A -A) 



X4 sin #d(9 
x 3 A x 



G 2 



tt/2 
J 2 



yiy2de 

2/3 Ay 



tt/2 

„ , n2 /" sin 2 6» cos 2 0d6> 
64^(^0-^) / - 





t/2 



(4.38c) 



21 A 2 



C 3 = 



V2yzde 
yiA y 



where we have defined 



X1X2X3X4, 



yi 2/2 2/32/4, 



2i2 2 2 3 2 4 , (4.39) 



and the factors Si, 2/i and Zi (i = 1, 2, 3, 4) are given by 
x\ = (A — /io) + (Ao— A) cos 2 6, £2 = (A— /i) + (Ao — A) cos 2 6>, 
x 3 = (A — fo) + (Ao — A) cos 2 6, X4 = (A— f) + (Ao — A) cos 2 9, 

J/1 = (/!— ^o) + (MO— H) COS 2 #, 2/2 = (/i — f ) + (/!() — /i) COS 2 0, 

2/ 3 = (/i — A ) + (/io — m) cos2 0> y4. = (p-X) + (po—[i) cos 2 0, 

2i = (l/— A()) + (f() — ^) cos 2 22 = (y— X)+(yo — u) cos 2 0, 

23 = {y — /Xq) + (vq — f) cos 2 0, 24 = (f— jt/) + (fo — cos 2 0. 

(4.40) 

For each i these factors follow from each other by cyclic 
permutation ofA^/i^z^^A and at the same time 
Ao — > Ho — ► v o — ► Ao- Half of the factors - all Xi, 2/1 and 
2/2 - are always positive, whereas the other factors are al- 
ways negative. The latter implies that one has to be careful 
with the signs of the factors under the square root when 
evaluating the single quadratures numerically. 

4-5.3 The complete elliptic integrals 

The two-dimensional homogeneous solutions F, G, H and I 
are given in l|4,90 in terms of the Legendre complete elliptic 
integrals E(m) and E'(m) = [E(m) — K(m)]/2m. Numer- 
ical routines for E(m) and K(m) (e.g., Press et al. 1999) 
generally require the argument to be < m < 1. In the 
allowed range of the confocal ellipsoidal coordinates, the ar- 
guments u 14. im and w H3.lt) I I become larger than unity. 
In these cases we can use transformations to express E(m) 
and K(m) in terms of E(l/m) and K(l/m) (e.g., Byrd & 
Friedman 1971). 

We prefer, however, to write the complete elliptic inte- 
grals as single quadratures similar to the above expressions 



for the hyperelliptic integrals. These quadratures can easily 
be evaluated numerically and apply to the full range of the 
confocal ellipsoidal coordinates. The resulting expressions 
for the two-dimensional homogeneous solutions are 



F(X,fx; X ,ho) 
G(X, (jl; A ,a*o) 
H(y, v; Ho, v ) 
v; Ho, fo) 



tt/2 

A— fi f Xid9 



TTVXo—fJ.0 J X2y/X\X2 ' 





A— h 



TT/2 



7r y Ao —/io 



sin 2 9d9 



2/V2/32/4 ' 



t/2 



H— v 



ty y Mo — vo 



(mo-m) 



2/2V2/12/2 ' 




(4.41) 



Again we have added two arguments to make the position 
of the unit source explicitly. We note that the homogeneous 
solutions A(X, /i; Ao, ho) and B(X, h\ Ao, po) for the disc case 
13. 4H are equivalent to F and G respectively. 



4.6 General triaxial solution 

We now construct the solution of the Jeans equations for 
triaxial Stackel models 14.21 . by superposition of singular 
solutions, which involve the homogeneous solution derived 
in the above. We match the solution to the boundary con- 
ditions at fj, — —a and v = —f3, and check for convergence 
of the solution when A — > 00. Next, we consider alternative 
boundary conditions and present the triaxial solution for a 
general finite region. We also show that the general solution 
yields the correct result in the case of thin tube orbits and 
the triaxial Abel models of Dejonghe & Laurent (1991). Fi- 
nally, we describe a numerical test of the triaxial solution to 
a polytrope model. 



4.6.1 Superposition of singular solutions 

Substitution of the functions A, B, C l|4.35[l and the func- 
tions F, G, H, I 14. All in expression 14. 61 . provides the 
three singular solutions of the system of simplified Jeans 
equations, with the right-hand side given by (|4.4p . We de- 
note these by S^ T (j = A, (i, v). The singular solutions of the 
two similar simplified systems, with the triplet of delta func- 
tions at the right-hand side of the first and third equation, 
Sl T and SJ 1 " then follow from 5*2 T by cyclic permutation. 
This gives 

51 A = B(u, A, fj,; vo, Ao, Ho)+G{v, A; u , Xo)S(ho-h) 

+H(X,(i; Ao, lM))5(vo — v) — 5(fM> — h)S(v — v), 
Si" = C(v, A, fx; vo, Ao, fM>)+I(X, fi; Ao, Ho)S(v — v) 
S" v = A{v, X, m vo, X , ho) + F(u, A; u , Xo)5(fJ,o— /i), (4.42a) 

52 X = A(X, h, v; A , fio, v )+F(X, /i; Ao, Ho)S{vo— v), 

= B(X, /Li, v; Ao, Mo, fo)+G(A, (i; A , ho)S{vq—v ) 
+H(fi, v\ fj,o, vo)S(Xo—X)—d(v —v)5(Xo — X), 
= C(A, h, v; A , Ho, vo)+I{h, v\ ho, v )S(X -X) (4.42b) 
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S3 = C(/i, v, A; no, iso, Xo)+I(v, X; i> , X )8(no-n), 
= A(n, v, A; fi , i/ , X )+F(n, Mo, f )<5(A -A) 
S^" = v, A; /io, fo, Ao) + G(/i, v; no, fo)5(Xo — X) 

+H(y, A; vo, Ao)5(/io-A«)-5(Ao-A))5(/io-/i). (4.42c) 

These singular solutions describe the contribution of a 
source point in (Ao,Mo, fo) to (A,/i, To find the solution 
of the full equations 1 1.21 , we multiply the singular solutions 
14.42all . Il4.42bll and 14.42cl by ffi(A , Mo, Vo), g2{Xo,po,vo) 
and g3(Ao, /-to, vo), respectively, so that the contribution from 
the source point naturally depends on the local density and 
potential (cf. eq. 14.31 1. Then, for each coordinate r = A, n, v, 
we add the three weighted singular solutions, and integrate 
over the volume fi, defined as 

fi = {(A , Mo, Mo) : A < A <oo, /i < Mo<-Q, v < v <-P}, (4.43) 

which is the three-dimensional extension of the integration 
domain D in Fig. [I] The resulting solution solves the in- 
homogeneous Jeans equations 11.21 . but does not give the 
correct values at the boundaries n = —a and v — —/3. They 
are found by multiplying the singular solutions 14.42 bt eval- 
uated at no = —a, and, similarly, the singular solutions 
14.42cl evaluated at Vq = —f3, by — iS MM (Ao, — a, vo) and 
— S vv (Xo, no, — P), respectively, and integrating in £1 over the 
coordinates that are not fixed. One can verify that this pro- 
cedure represents the boundary values correctly. The final 
result for the general solution of the Jeans equations 114.211 
for triaxial Stackel models is 

oo — a —f3 ^ 

S TT (X,n,v) =jdXoJdfj,oJdvo^>2 gi(Xo,no,vo)Sl T (X,n,v;X ,no,vo) 

A m v 8=1 
-0 oo 

-JduojdXo Sfinfioj—apo) S2 T (X,n,u;X ,—a,vo) 

v A 
oo —a 

-JdXojdno S vv (Xo,no,-P) 5 , 3 rT (A,/i,^;A ,Mo,-/3),(4.44) 

A j.L 

where r = (A, n, v \ This gives the stresses everywhere, once 
we have specified S^iX, — a, v) and S VV (X, n, — P). At both 
boundaries /i = —a and v — —/3, the three stress compo- 
nents are related by a set of two Jeans equations, i.e., (14. 211 
evaluated at p = ~ a an d v = — y3 respectively. From 
we know that the solution of these two-dimensional systems 
both will involve a (boundary) function of one variable. We 
need this latter freedom to satisfy the continuity conditions 
12.171 1. This means it is sufficient to specify any of the three 
stress components at fj, = —a and v — —p. 

4-6.2 Convergence of the general triaxial solution 

As in § H3.1.4l|3.2.7l and 13.41 we suppose G(t) = 0{t 6 ) when 
t — *• oo, with 8 in the range [—5,0). This implies that the 
potential Vs EHJ is also 0(t s ). We assume that the density 
p, which does not need to be the density ps which generates 
Vs, is of the form N(ji, u)X~ s ^ 2 when A — * 00. In the special 
case where p — ps, we have s < 4 except possibly along the 
2-axis. When s = 4 the models remain flattened out to the 
largest radii, but when s < 4 the function N(p, v) 1 in 
the limit A — *■ 00 (de Zeeuw et al. 1986). 

From the definition (14.311 . we find that gi(Xo, no, vo) = 



£>(Xo~ ) as A co, while g2(X , no, Vo) and g s (X ,no, vo) 
are larger and both O(A s ^ 2 ). To investigate the behaviour 
of the singular solutions (14.421 at large distance, we have 
to carefully analyse the complete hyperelliptic 14.3811 and 
elliptic 14.411 integrals as Ao — > 00. This is simplified by 
writing them as Carlson's 7?- functions (Carlson 1977). We 
finally find for the singular solutions that S{ T = 0(1) when 
Ao — > 00, whereas SJ 1 ^ and SJ T are smaller and 0(Xq 1 ), with 
t — A, n, v. This shows that for the volume integral in the 
triaxial solution 14.441 to converge, we must have S — s/2 + 
1 < 0. This is equivalent to the requirement s > 28 + 2 we 
obtained in H3.4l for the limiting cases of prolate and oblate 
potentials and for the large radii limit with scale-free DF. 
From the convergence of the remaining two double integrals 
in l|4.44|l . we find that the boundary stresses S , )JM (A, — a, v) 
and S VV (X, n, — j9) cannot exceed O(l) when A — > 00. 

The latter is in agreement with the large A behaviour of 
S TT {X, n, v ) that follows from the volume integral. The sin- 
gular solutions Si X = 0(1) (i = 1, 2, 3) when A — > 00, larger 
than and S"" , which are all C(A _1 ). Evaluating the 
volume integral at large distance then gives S TT (X, n, v) = 
0{X s ~ a / 2+1 ), i.e., not exceeding 0(1) if the requirement 
s > 28 + 2 is satisfied. We obtain the same behaviour and 
requirement from the energy equation 12.101 . 

We conclude that for the general triaxial well 
as for the limiting cases with a three-dimensional shape, 
the stress components T TT (X,n, v) are 0(X S ~ S ^ 2 ) at large 
distance, with the requirement that s > 28 + 2 for — | < 
8 < 0. We obtained the same result for the stresses in the 
disc case, except that then s > 28 + 1. Both the three- 
dimensional and two-dimensional requirements are met for 
many density distributions p and potentials Vs of interest. 
They do not break down until the isothermal limit 5^0, 
with s = 1 (disc) and s = 2 (three-dimensional) is reached. 



4-6.3 Alternative boundary conditions 

Our solution for the stress components at each point (A, p, v) 
in a triaxial model with a Stackel potential consists of the 
weighted contribution of all sources outwards of this point. 
Accordingly, we have integrated with respect to Ao, no an d 
vo, with lower limits the coordinates of the chosen point 
and upper limits 00, —a and —j3, respectively. To obtain 
the correct expressions at the outer boundaries, the stresses 
must vanish when A — > 00 and they have to be specified at 
n = —a and v — —j3. 

The integration limits A, n an d v are fixed, but for the 
other three limits we can, in principle, equally well choose 
—a, — P and —7 respectively. The latter choices also imply 
the specification of the stress components at these bound- 
aries instead. Each of the eight possible combinations of 
these limits corresponds to one of the octants into which the 
physical region —7 < v < — P < Mo < ~ a < Ao < 00 is split 
by the lines through the point (A, n, v). By arguments simi- 
lar to those given in H3.3I one may show that in all octants 
the expressions 14.351 for A, B, C, and (PHI for F, G, H, I 
are equivalent. Hence, again the only differences in the sin- 
gular solutions are due to possible changes in the sign of the 
step-functions, but the changes in the integration limits can- 
cel the sign differences between the corresponding singular 
solutions. However, as in H3.3l for the two-dimensional case, 
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it is not difficult to show that while switching the boundary 
conditions fi and v is indeed straightforward, the switch from 
A — > oo to A = —a again leads to solutions which generally 
have the incorrect radial fall-off, and hence are non-physical. 



4.6.4 Triaxial solution for a general finite region 

If we denote non-fixed integration limits by A e , fi e and v e 
respectively, we can write the triaxial solution for a general 
finite region as 

A e Me V & g 

A H v 1=1 

Me f e 

d/iajduo S\\(X e ,(j,o,vo) Sl T (X,fi,u;X e ,no,i'o) 

M V 
V e Ae 

d^o^dAo S^XcHefVo) SI t {X,h,v\Xo,Hz,vo) 

A 

e Me 

dAo dfj,o Sw(Xo,no,Ve) Sj T (A,/n,^;Ao,/io,^e), (4.45) 



v A 
A e M 



with, as usual, r = A, fx, v. The weight functions gi (i — 
1, 2, 3) are defined in 14.31 and the singular solutions SJ T are 
given by (14.4211 . The non- fixed integration limits are chosen 
in the corresponding physical ranges, i.e., X e 6 [—a, 00], 
He £ [— ft, —a] and v e g [—7, — ft], but A e / —a (see H4.6.31 . 
The solution requires the specification of the stress compo- 
nents on the boundary surfaces X = X e , fi = fi e and v = z/ e . 
On each of these surfaces the three stress components are re- 
lated by two of the three Jeans equations 14.21 and the conti- 
nuity conditions 112. 1711 . Hence, once one of the stress compo- 
nents is prescribed on three boundary surfaces, the solution 
14.441 1 yields all three stresses everywhere in the triaxial 
Stackel galaxy. The stresses on the remaining three bound- 
ary surfaces then follow as the limits of the latter solution. 



DF is now fixed, and can be found by solving a system of 
linear equations, starting from the outside (A — > 00). 

The total DF is the sum of the DFs of the four orbit 
families, and is hence highly non-unique. All these DFs give 
rise to a range of stresses Taa, ?mm> T vv , and our solution of 
the Jeans equations must be sufficiently general to contain 
them as a subset. This is indeed the case, as we are allowed 
to choose the stresses on the special surfaces v — —ft and 
H — —a. However, not all choices will correspond to physical 
DFs. The requirement T TT > is necessary but not sufficient 
for the associated DF to be non-negative everywhere. 



4.6.6 The general solution for thin tube orbits 

For each of the three tube families in case of infinitesimally 
thin orbits one of the three stress components vanishes ev- 
erywhere (see H2.5.61 . We are left with two non-zero stress 
components related to the density and potential by three 
reduced Jeans equations 14.21 . We thus have subsidiary con- 
ditions on the three right hand side terms gi, g2 and 173. 

HZ92 solved for the two non-trivial stresses and showed 
that they can be found by single quadratures (with inte- 
grands involving no worse than complete elliptic integrals), 
once the corresponding stress had been chosen at v = —ft 
(for I- and O-tubes) or at /1 = —a (for S-tubes). 

By analogy with the reasoning for the thin tube orbits in 
the disc case ( H3.4.41 . we can show that for each of the three 
tube families in the case of thin orbits the general triaxial so- 
lution 14.451 gives the stress components correctly. Consider, 
e.g., the thin I-tubes, for which S MM = 0. Apply the latter 
to 14. 4511 . substitute for gi, g^ and 513 the subsidiary condi- 
tions that follow from the reduced Jeans equations 14.21 and 
substitute for the singular solutions the expressions 14.421 . 
After several partial integrations, we use that the homoge- 
neous solutions A, B and C solve a homogeneous system 
similar to 14.111 . but now with respect to the source point 
coordinates (Ao,po,^o) 



4-6.5 Physical solutions 

Statler (1987) and HZ92 showed that many different DFs 
are consistent with a triaxial density p in the potential Vs- 
Specifically, the boundary plane v — —ft, i.e., the area out- 
side the focal hyperbola in the (x, z)-plane (Fig. is only 
reached by inner (I) and outer (O) long-axis tube orbits. A 
split between the contribution of both orbit families to the 
density in this plane has to be chosen, upon which the DF 
for both the I and O orbits is fixed in case only thin tubes 
are populated, but many other possibilities exist when the 
full set of I- and O-orbits is included. For each of these DFs, 
the density provided by the I- and O-tubes can then in prin- 
ciple be found throughout configuration space. In the area 
outside the focal ellipse in the (y, z)-plane (/i = —a), only 
the O-tubes and S-tubes contribute to the density. Subtract- 
ing the known density of the O-orbits leaves the density to 
be provided by the S-tubes in this plane, from which their 
DF can be determined. This is again unique when only thin 
orbits are used, but is non-unique otherwise. The density 
that remains after subtracting the I-, O-, and S-tube densi- 
ties from p must be provided by the box (B) orbits. Their 



dB(v,X,n;v ,X ,no) _ A(X,n,u;X ,no,ni) C{pL,v,X;p ,v ,Xa) 



dX 



2(A — Ho) 



2(Xo-u ) 



(4.46) 



where other relations follow by cyclic permutation of A — > 
H — > v — > A and Ao — ► Ho — * v o — > -^o- And similar for the 
two-dimensional homogeneous solutions F, G, H and / the 
relations follow from 



<9G(/i, A; ho, Ao) _ F(X,h;X ,ho) 



dXo 

8H(h, v\ho,vq) 
dpo 



2(Ao — ho) 

I(v,H; vq,Ho) 
2(ho — fo) 



(4.47) 



It indeed turns out that for 5^ (A, v) all terms cancel on 
the right hand side of 14.451 . The terms that are left in 
the case of Saa and S vv are just eq. I)4.2a|l integrated with 
respect to A and eq. 14.2cl integrated with respect to v, 
respectively, and using that S^^ = 0. A similar analysis as 
above shows that also for thin O- and S-tubes - for which 
S\\ = in both cases - the general triaxial solution yields 
the correct result. 
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4-6.7 Triaxial Abel models 

For a galaxy with a triaxial potential of Stackel form, the DF 
is a function of the three exact isolating integrals of motion, 
/(x, v) = f(E,l2,l3) (see also >I2.21 . The expressions for 
E, I2 and I3 in terms of the phase-space coordinates (x, v) 
can be found in e.g. Z85. We can thus write the velocity 
moments of the DF as a triple integral over E, I2 and I3. 
Assuming that the DF is function of only one variable 



S = E + wh + uh 



(4.48) 



with w and u constants, Dejonghe & Laurent (1991) show 
that the triple integration simplifies to a one-dimensional 
Abel integration over S. Even though a DF of this form can 
only describe a self-consistent model in the spherical case 
(ellipsoidal hypothesis, see, e.g., Eddington 1915), the Jeans 
equations do not require self-consistency. 

The special Abel form results in a simple analytical re- 
lation between the three stress components (Dejonghe & 
Laurent 1991, their eq. [5.6]) 



Tfiy. — Twa^v/axv, 
with 



Ti/v — TxxOj flu I ' ct^x, 



(4.49) 



aw = {l — a) + (a + a)(r + a)w — (a+-/)(r+j)u, (4.50) 
and cr, t = A, p, v. With these relations we find that 

(4.51) 



Txx~T u 



Txx daxv 



-T vv Txx daxu 



X — p 



ax v 



dX ' 



X-u 



axu 



dX 



The first Jeans equation l|2.16a[l now becomes a first-order 
partial differential equation for Txx- This equation can be 
solved in a straightforward way and provides an elegant and 
simple expression for the radial stress component 



axg^axov OVs 



(4.52) 



T X x{X,p,v) = ax "V ax ev T\x{X e ,(J,,v) 
V a-x^a-Xv 

x c 

+ JdX 

x 

The expressions for T MM and T vv follow by application of the 
ratios 14. 491 . If we let the boundary value A e — > 00, the first 
term on the right-hand side of 14.521 vanishes. 

The density p, which does not need to be the density ps 
which generates Vs, is of the Abel form as given in eq. (3.11) 
of Dejonghe & Laurent (1991). If we substitute this form in 
14. 521 . we obtain, after changing the order of integration and 
evaluating the integral with respect to A, again a single Abel 
integral that is equivalent to the expression for Txx that 
follows from eq. (3.10) of by Dejonghe & Laurent (1991). 
Using the relations 14.491 and the corresponding subsidiary 
conditions for gi, Q2 and gs, it can be shown that the general 
triaxial solution 14.451 gives the stress components correctly. 

4-6.8 Numerical test 

We have numerically implemented the general triaxial solu- 
tion 14.451 . and tested it on a polytrope dynamical model, 
for which the DF depends only on energy E as f(E) oc 
_E" _3//2 , with n > i. Integration of this DF over velocity v, 
with E = — V — \v 2 for a potential V < 0, shows that the 
density p oc (-V) n (e.g., Binney & Tremaine 1987, p. 223). 
This density is not consistent with the Stackel potentials we 



use but, as noted in H2.3I the Jeans equations do not require 
self-consistency. The first velocity moments and the mixed 
second moments of the DF are all zero. The remaining three 
moments all equal —V/(n+ 1), so that the isotropic stress 
of the polytrope model T pol oc (— V) n+1 . 

We take the potential V to be of Stackel form Vs 12.31 . 
and consider two different choices for G(r) in 12.41 . The first 
is the simple form G(t) = — GAf/ (\/t"+\/~ ct) that is related 
to Henon's isochrone (de Zeeuw & Pfenniger 1988). The sec- 
ond is the form for the perfect ellipsoid, for which G(t) is 
given in Z85 in terms of complete elliptic integrals. The par- 
tial derivatives of Vs(X, p, v), that appear in the weights <?i, 
g2 and <?3, can be obtained in terms of G(t) and its derivative 
in a straightforward way by using the expressions derived by 
de Zeeuw et al. (1986). 

The calculation of the stresses is done in the following 
way. We choose the polytrope index n, and fix the triaxial 
Stackel model by choosing a, /3 and 7. This gives T vo \. Next, 
we obtain successively the stresses Txx, T^y, and T vv from 
the general triaxial solution 14.451 by numerical integration, 
where the relation between S TT and T TT is given by 14.11 . We 
first fix the upper integration limits A e , p e and v e . All inte- 
grands contain the singular solutions 14.421 . that involve the 
homogeneous solutions A, B, G, F, G, H and 7, for which 
we numerically evaluate the single quadratures (eq. 14.351 . 
14.381 and 14.411 1. The weights gi, g 2 and g 3 14.31 involve 
the polytrope density and Stackel potential. This leaves the 
boundary stresses in the integrands, for which we use the 
polytrope stress T po i that follows from the choice of the DF, 
evaluated at the corresponding boundary surfaces. We then 
evaluate the general solution away from these boundaries, 
and compare it with the known result. 

We carried out the numerical calculations for different 
choices of n, a, j3 and 7 and at different field points (A, [i, v). 
In each case the resulting stresses Txx, and T„„ - inde- 
pendently calculated - were equivalent to high precision and 
equal to T po \. This agreement provides a check on the accu- 
racy of both our formulae and their numerical implementa- 
tion, and demonstrates the feasibility of using our methods 
for computing triaxial stress distributions. That will be the 
subject of a follow-up paper. 



5 DISCUSSION AND CONCLUSIONS 

Eddington (1915) showed that the velocity ellipsoid in a tri- 
axial galaxy with a separable potential of Stackel form is 
everywhere aligned with the confocal ellipsoidal coordinate 
system in which the equations of motion separate. Lynden- 
Bell (1960) derived the three Jeans equations which relate 
the three principal stresses to the potential and the density. 
They constitute a highly-symmetric set of first-order par- 
tial differential equations in the three confocal coordinates. 
Solutions were found for the various two-dimensional limit- 
ing cases, but with methods that do not carry over to the 
general case, which, as a consequence, remained unsolved. 

In this paper, we have introduced an alternative solu- 
tion method, using superposition of singular solutions. We 
have shown that this approach not only provides an elegant 
alternative to the standard Riemann-Green method for the 
two-dimensional limits, but also, unlike the standard meth- 
ods, can be generalised to solve the three-dimensional sys- 
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tem. The resulting solutions contain complete (hyper)elliptic 
integrals which can be evaluated in a straightforward way. 
In the derivation, we have recovered (and in some cases cor- 
rected) all previously known solutions for the various two- 
dimensional limiting cases with more symmetry, as well as 
the two special solutions known for the general case, and 
have also clarified the restrictions on the boundary values. 
We have numerically tested our solution on a polytrope 
model. 

The general Jeans solution is not unique, but requires 
specification of principal stresses at certain boundary sur- 
faces, given a separable triaxial potential, and a triaxial den- 
sity distribution (not necessarily the one that generates the 
potential) . We have shown that these boundary surfaces can 
be taken to be the plane containing the long and the short 
axis of the galaxy, and, more specifically, the part that is 
crossed by all three families of tube orbits and the box or- 
bits. This is not unexpected, as HZ92 demonstrated that the 
phase-space distribution functions of these triaxial systems 
are defined by specifying the population of each of the three 
tube orbit families in a principal plane. Once the tube orbit 
populations have been defined in this way, the population 
of the box orbits is fixed, as it must reproduce the density 
not contributed by the tubes, and there is only one way to 
do this. While HZ92 chose to define the population of inner 
and outer long axis tubes in a part of the {x, z)-plane, and 
the short axis tubes in a part of the (y, z)-plane, it is in fact 
also possible to specify all three of them in the appropriate 
parts of the (x, z)-plane, just as is needed for the stresses. 

The set of all Jeans solutions 14.451 contains all the 
stresses that are associated with the physical distribution 
functions / > 0, but, as in the case of spherical and axisym- 
metric models, undoubtedly also contains solutions which 
are unphysical, e.g., those associated with distribution func- 
tions that are negative in some parts of phase space. The 
many examples of the use of spherical and axisymmetric 
Jeans models in the literature suggest nevertheless that the 
Jeans solutions can be of significant use. 

While triaxial models with a separable potential do not 
provide an adequate description of the nuclei of galaxies with 
cusped luminosity profiles and a massive central black hole, 
they do catch much of the orbital structure at larger radii, 
and in some cases even provide a good approximation of 
the galaxy potential. The solutions for the mean streaming 
motions, i.e., the first velocity moments of the distribution 
function, are quite helpful in understanding the variety of 
observed velocity fields in giant elliptical galaxies and con- 
straining their intrinsic shapes (e.g., Statler 1991, 1994b; 
Arnold et al.1994; Statler, DeJonghe & Smecker-Hane 1999; 
Statler 2001). We expect that the projected velocity disper- 
sion fields that can be derived from our Jeans solutions will 
be similarly useful, and, in particular, that they can be used 
to establish which combinations of viewing directions and 
intrinsic axis ratios are firmly ruled out by the observations. 
As some of the projected properties of the Stackel models 
can be evaluated by analytic means (Franx 1988) , it is possi- 
ble that this holds even for the intrinsic moments considered 
here. Work along these lines is in progress. 

The solutions presented here constitute a significant 
step towards completing the analytic description of the prop- 
erties of the separable triaxial models, whose history by now 
spans more than a century. It is remarkable that the entire 



Jeans solution can be written down by means of classical 
methods. This suggests that similar solutions can be found 
for the higher dimensional analogues of 12.161 . most likely 
involving hyperelliptic integrals of higher order. It is also 
likely that the higher-order velocity moments for the sepa- 
rable triaxial models can be found by similar analytic means, 
but the effort required may become prohibitive. 



ACKNOWLEDGMENTS 

This paper owes much to Donald Lynden-BelPs enthusiasm 
and inspiration. This work begun during a sabbatical visit 
by CH to Leiden Observatory in 1992, supported in part by a 
Bezoekersbeurs from NWO, and also by NSF through grant 
DMS 9001404. CH is currently supported by NSF through 
grant DMS 0104751. This research was supported in part 
by the Netherlands Research School for Astronomy NOVA. 
The authors gratefully acknowledge stimulating discussions 
with Wyn Evans during the initial phases of this work. 



REFERENCES 

Abramowitz M., Stegun I. A., 1965, Handbook of Mathematical 

Functions. New York, Dover Publications 
Arnold R., 1995, MNRAS, 276, 293 

Arnold R., dc Zeeuw P. T., Hunter C, 1994, MNRAS, 271, 924 

Bacon R., Simien F., Monnet G., 1983, A&A, 128, 405 

Bak J., Statler T. S., 2000, AJ, 120, 110 

Binney J., 1976, MNRAS, 177, 19 

Binney J., 1978, MNRAS, 183, 501 

Binney J., 1982, MNRAS, 201, 15 

Binney J., Tremaine S., 1987, Galactic Dynamics. Princeton, NJ, 

Princeton University Press 
Bishop J. L., 1986, ApJ, 305, 14 
Bishop J. L., 1987, ApJ, 322, 618 

Byrd P. F., Friedman M. D., 1971, Handbook of Elliptic Integrals 
for Engineers and Scientists. Springer- Verlag, Berlin Heidel- 
berg New York, 2nd revised ed. 

Carlson B. C, 1977, Special Functions of Applied Mathematics. 
Academic Press, New York San Francisco London 

Chandrasekhar S., 1939, ApJ, 90, 1 

Chandrasekhar S., 1940, ApJ, 92, 441 

Conway J., 1973, Functions of one complex variable. New York, 
Springer- Verlag 

Copson E. T., 1975, Partial Differential Equations. Cambridge, 

Cambridge University Press 
de Zeeuw P. T., 1985a, MNRAS, 216, 273 [Z85] 
de Zeeuw P. T., 1985b, MNRAS, 216, 599 
de Zeeuw P. T., Hunter C, 1990, ApJ, 356, 365 
de Zeeuw P. T., Hunter C, Schwarzschild M., 1987, ApJ, 317, 

607 

de Zeeuw P. T., Peletier R., Franx M., 1986, MNRAS, 221, 1001 
de Zeeuw P. T., Pfenniger D., 1988, MNRAS, 235, 949 
Dejonghe H., de Zeeuw P. T., 1988, ApJ, 333, 90 
Dejonghe H., Laurent D., 1991, MNRAS, 252, 606 
Eddington A. S., 1915, MNRAS, 76, 37 
Evans N., 1990, Intern. J. Computer Math., 34, 105 
Evans N. W., Carollo C. M., de Zeeuw P. T., 2000, MNRAS, 318, 
1131 

Evans N. W., de Zeeuw P. T., 1992, MNRAS, 257, 152 
Evans N. W., Lynden-Bell D., 1989, MNRAS, 236, 801 [EL89] 
Evans N. W., Lynden-Bell D., 1991, MNRAS, 251, 213 
Franx M., 1988, MNRAS, 231, 285 



General solution Jeans equations 27 



Gebhardt K., Richstone D., Tremaine S., Lauer T. D., Bender R., 
Bower G., Dressier A., Faber S., Filippenko A. V., Green R., 
Grillmair C, Ho L. C, Kormendy J., Magorrian J., Pinkncy 
J., 2003, ApJ, in press 

Gerhard O. E., 1993, MNRAS, 265, 213 

Goldstein H., 1980, Classical Mechanics. London, Addison- Wesley 
Hunter C, de Zeeuw P. T., 1992, ApJ, 389, 79 [HZ92] 
Hunter C, de Zeeuw P. T., Park C, Schwarzschild M., 1990, ApJ, 
363, 367 

Jeans J. H., 1915, MNRAS, 76, 70 

Kuzmin G., 1973, in Proc. Ail-Union Conf., Dynamics of Galaxies 
and Clusters, ed. T.B. Omarov (Alma Ata: Akad. Nauk Kaza- 
khskoj SSR), 71 (English transl. in IAU Symp. 127, Structure 
and Dynamics of Elliptical Galaxies, ed. P.T. de Zeeuw [Dor- 
drecht: Reidel], 553) 

Lyndcn-Bell D., 1960, PhD thesis, Cambridge University 

Lynden-Bell D., 1962a, MNRAS, 124, 1 

Lynden-Bell D., 1962b, MNRAS, 124, 95 

Mathieu, A., Dejonghe, H., 1999, MNRAS, 303, 455 

Oberhettinger F., Badii L., 1973, Tables of Laplace Transforms. 
New York, Springer- Verlag 

Press W. H., Teukolsky S. A., Vettering W. T., Flannery B. P., 
1992, Numerical Recipes. Cambridge Univ. Press, Cambridge 

Qian E. E., de Zeeuw P. T., van der Marel R. P., Hunter C, 1995, 
MNRAS, 274, 602 

Rix H., de Zeeuw P. T., Cretton N., van der Marel R. P., Carollo 
C. M., 1997, ApJ, 488, 702 

Robijn F. H. A., de Zeeuw P. T., 1996, MNRAS, 279, 673 

Schwarzschild M., 1979, ApJ, 232, 236 

Schwarzschild M., 1993, ApJ, 409, 563 

Stackel P., 1890, Math. Ann., 35, 91 

Stackel P., 1891, Uber die Integration der Hamilton- Jacobischen 
Differential gleichung mittclst Separation der Variabeln. Ha- 
bilitationsschrift, Halle 

Statler T. S., 1987, ApJ, 321, 113 

Statler T. S., 1991, AJ, 102, 882 

Statler T. S., 1994a, ApJ, 425, 458 

Statler T. S., 1994b, ApJ, 425, 500 

Statler T. S., 2001, AJ, 121, 244 

Statler T. S., Fry A. M., 1994, ApJ, 425, 481 

Statler, T. S., Dejonghe, H., Smecker-Hane, T., 1999, AJ, 117, 
126 

Strauss W. A., 1992, Partial Differential Equations. New York, 
John Wiley 

van der Marel R. P., Cretton N., de Zeeuw P. T., Rix H., 1998, 
ApJ, 493, 613 

Verolme E. K., Cappellari M., Copin Y., van der Marel R. P., 
Bacon R., Bureau M., Davies R., Miller B., de Zeeuw P. T., 
2002, MNRAS, 335, 517 

Weinacht J., 1924, Math. Ann., 91, 279 

Zhao H. S., 1996, MNRAS, 283, 149 



APPENDIX A: SOLVING FOR THE 
DIFFERENCE IN STRESS 

We compare our solution for the stress components T\\ and 
T/jm with the result derived by EL89. They combine the two 
Jeans equations (12.2511 into the single equation 
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for the difference A = T\\ — T MM of the two stress compo- 
nents. Eq. 1A11 is of the form 
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Figure Al. The physically relevant region of the (A, p)-plane for 
the determination of the Riemann— Green function G, overlayed 
with the new coordinates £ and n IA5I . The dot marks the source 
point of the Riemann-Green function G at (Ao, po). This function 
is non-zero only in the shaded region, which denotes the domain 
of influence in the (A, p)-plane of that source point. Fie. 131 on the 
other hand shows the (Ao, po)- plane. It is relevant to the solution 
for the stress at a single held point (A, p). The hatched region D 
of Fie. Bl shows the domain of dependence of the held point, that 
is the portion of the source plane on which the solution at the 
field point depends. 



where C is the adjoint operator defined in eq. (13.61 . As in 
JO eq. ED can be solved via a Riemann-Green function. 



Al The Green's function 

In order to obtain the Riemann-Green function Q* for the 
adjoint operator C* , we use the reciprocity relation (Copson 
1975, §5.2) to relate it to the Riemann-Green function Q, de- 



rived in fln~21 for C. With ci = c 2 

Q*(X,p;X ,po) = G{X ,po;X,p) 

1 

/ Xp-pLQ 

~ V x—p 



2 in this case, we get 



Fi(-i,f;l;to), 



(A3) 



where w as defined in i3~ToTl . EL89 seek to solve eq. llAH 
using a Green's function G which satisfies the equation 



C*G = 5{X -X)S(p -p). 



(A4) 



£*A 



(A2) 



That they impose the same boundary conditions that we do 
is evident from their remark that, if C* were the simpler op- 
erator d 2 /dXdp, G would be TL{Xo—X)7i(po—p). This is the 
same result as would be obtained by the singular solution 
method of H3.2I which, as we showed there, is equivalent to 
the Riemann-Green analysis. Hence their G should match 
the Q* of eq. 1A3II . We show in flA3l that it does not. 
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A2 Laplace transform 

We use a Laplace transform to solve (IA4I because the re- 
quired solution is that to an initial value problem to which 
Laplace transforms are naturally suited. The PDE is hyper- 
bolic with the lines A = const and [i = const as charac- 
teristics, and its solution is non-zero only in the rectangle 
bounded by the characteristics A = Ao and /i = /Jo, and the 
physical boundaries A = —a and fi = — /3 (Fig. lAH . We 
introduce new coordinates 

£ = (A-a*)/V2, V = -(A+m)A/2, (A5) 
so that eq. 1A4I simplifies to 

d 2 G 8 2 G d (G\ nsfi . . . , . . r ., 

9^- W -9i{j)=^-^(v-Vo), (A6) 

where £o = (-^o — A t o)/\/2 and 770 = — (Ao + /tto)/V2 are the 
coordinates of the source point. The factor of 2 arises from 
the transformation of the derivatives; the product of the 
delta functions in l!A4(l transforms into that of l|A6(l because 
the Jacobian of the transformation 1A5H is unity. The reason 
for our choice of r\ is that G = for r\ < 770, that is A + 
H > Ao + /io • Hence rj is a time- like variable which increases 
in the direction in which the non-zero part of the solution 
propagates. We take a Laplace transform in fj — r\ — r/o, and 
transform G(£, rf) to 



C*G : 



G(£,p) 



e-^G(^,fj)dv- 



(A7) 



There are two equally valid ways of taking proper account 
of the 5(rj — r/o) in taking the Laplace transform of equation 
dXet . One can either treat it as 5(77—0+) , in which case it has 
a Laplace transform of 1, or one can treat it as 8(fj — 0— ), 
in which case it contributes a unit initial value to dG/drj 
which must be included in the Laplace transform of d 2 G/drf 
(Strauss 1992). Either way leads to a transformed equation 
for G(£,p) of 



2 A 
p G - 



d 2 G 
d£2" 



m - Co) 



(A8) 



The homogeneous part of eq. <A8t is the modified Bessel 
equation of order one in the variable p£. Two independent 
solutions are the modified Bessel functions 7i and Ki. The 
former vanishes at £ = and the latter decays exponentially 
as £ — > oo. We need G to decay exponentially as £ — > oo 
because G(£, 77) vanishes for 77 < £ — £0, and hence its Laplace 
transform G is exponentially small for large f . We also need 
G to vanish at £ = where A = fi. The focus at which 
A = H = —a is the only physically relevant point at which 
£ = 0. It lies on a boundary of the solution region in the 
Ao — > —a limit (Fig. IA11 . The focus is a point at which 
the difference A between the stresses vanishes, and hence G 
and G should vanish there. The delta function in eq. 1A8I 
requires that G be continuous at £ = £0 and that dG/d£ 
decrease discontinuously by 2 as £ increases through £ = £o- 
Combining all these requirements, we obtain the result 

d(f > p) = ( 2&1Cl<pC)/l<P&) ' ^<°°> (A9) 

We use the Wronskian relation Ii(x)K[(x) — I[(x)Ki(x) = 
— 1/x (eq. [9.6.15] of Abramowitz & Stegun 1965) in cal- 
culating the prefactor of the products of modified Bessel 



functions. The inversion of this transform is obtained from 
formula (13.39) of Oberhettinger & Badii (1973) which gives 



|?o-C|<77<£o+£, 



(A10) 



(All) 



.0, -oo<?7<|£o-£|, 

we have (cf. eq. |3.16p 

w = f-i^-t) 2 = (Aq-AXmo-zj) 
4^oC (A -Mo)(A-/x)' 

The second case of eq. <A10t shows that G does indeed va- 
nish outside the shaded sector A < Ao , < Ho ■ The first case 
shows that it agrees with the adjoint Riemann-Green func- 
tion Q* of d A3I > which was derived from the analysis of H3.ll 

A3 Comparison with EL89 

EL89 use variables s = — 77 and t = £, whereas we avoided 
using t for the non-time-like variable. They consider the 
Fourier transform 



G{t,k) = 



-ikfj 



G(C,77)d77. 



(A12) 



Because G = for 77 < 0, we can rewrite our Laplace 
transform as their Fourier transform. Setting p = — ik 
gives G(£, k) = iG(£, — ik), and using the formulas Ii(x) = 
— Ji(ix) and Kx(x) = ^itiH^ (ix), eq. (IA9II yields 



G(£,fc) 



7rt£oi*r(fcf) Ji(fcfo), Co<C<oo, 



(A13) 



This formula differs from the solution for the Fourier 
transform given in eq. (70) of EL89. The major difference is 
that their solution has Hankel functions of the second kind 
H[ 2 \kt) = H[ 2) (k£) where ours has Ji Bessel functions. 
Consequently their solution has an unphysical singularity at 
t = £ = 0, and so, in our opinion, is incorrect. Our solution, 
which was devised to avoid that singularity, gives a result 
which matches that derived by Riemann's method in H3.ll 



A4 The solution for A 

The solution for A using the adjoint Riemann-Green func- 
tion is given by eq. 13.14H with Q replaced by Q* and the 
sign of C2 changed for the adjoint case (Copson 1975). The 
hypergeometric function of eq. (IA3I for Q* is expressible in 
terms of complete elliptical integrals as 
2 , 



2*1 (- 



1; w) = — [E(w) +2wE'(w)]. 



(A14) 



Hence, the solution for the difference A between the two 
principal stresses is given by 



A(X,tx) = 



7r(A — /J,) ' 



oo —a. 



dA 



E(w) + 2wE'( 



w 

M0 = -<* 



d\ dfi 9mo d\ J 
jL[(X +a)iA(X ,-aA. (A15) 
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The determined reader can verify, after some manipulation, 
that this expression is equivalent to the difference between 
the separate solutions l|M.21a|l and (|3. 21 bfl . derived in A3. II 

Note added in manuscript 

We agree with the amendment to our method of solution 
for A given in Appendix IA4I Our Green's function, while 
solving the differential equation, had the wrong boundary 
conditions. 

N.W. Evans & D. Lynden-Bell 
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